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Preface 


In  1978  an  article  in  the  Journal  of  Computational 
Physics  (Ref  6)  evaluated  numerical  methods  for  elliptical 
partial  differential  equations.  The  results  showed  the 
method  of  weighted  residuals  used  with  finite  element 
analysis  yielded  more  accurate  results  in  less  computer 
time  than  finite  difference  methods.  Dr.  Kaplan  of  the 
Air  Force  Institute  of  Technology,  Physics  Department,  and 
Dr.  Kessler  of  the  Air  Force  Materials  Laboratory  decided 
to  investigate  the  method  of  weighted  residuals  for  numer¬ 
ical  solution  of  the  heat  equation.  Specific  application 
would  eventually  be  to  the  transient  heat  equation  for 
bodies  of  axial  symmetry  with  the  goal  of  reduced  computer 
time  for  solution.  This  thesis  represents  a  first  step  in 
that  investigation  by  studying  the  method  of  weighted 
residuals  applied  on  the  whole  domain  of  interest  as  a 
numerical  method. 

Thanks  are  due  to  Dr.  Kaplan  for  his  unfailing  encour¬ 
agement  and  guidance  during  the  research  quarters.  Thanks 
are  also  due  to  Sally  Lindsay  who  typed  this  manuscript. 

Her  experience  and  professional  approach  did  much  to  present 
the  finished  product  in  its  best  light.  Last,  my  greatest 
thanks  and  appreciation  are  due  to  my  wife  and  infant 
daughter  for  their  personal  sacrifices  during  a  time  which 
proved  to  be  one  of  great  personal  adjustment  and  turmoil. 

Robert  E.  Naegeli 


ii 


Contents 

Page 

Preface .  ii 

List  of  Figures .  iv 

List  of  Tables .  v 

Abstract .  vi 

I.  Introduction  .  1 

II.  Theory .  3 

The  Method  of  Weighted  Residuals  .  3 

Types  of  Method  Application  .  6 

Types  of  Weighting  Criteria  .  7 

Trial  Functions .  9 

Transient  Problems  .  11 

Numerical  Solution  for  MWR .  17 

Steady  State  .  17 

Transient .  19 

The  Finite  Difference  Method  .  21 

III.  The  Problem  Set . 23 

Steady  State  Problem  .  23 

Transient  Problem  .  28 

IV.  Results  and  Comparisons .  36 

Measurement  of  Accuracy  .  37 

Steady  State  Problem  .  39 

Transient  Problem  .  47 

V.  Conclusions .  57 

Conclusions  and  Limitations  .  57 

Recommendations  .  60 

Bibliography  .  66 

Appendix  A:  Error  and  Solution  Time  Tables  .  68 

Appendix  B:  Computer  Programs  .  75 


iii 


Figure 

1 


List  of  Figures 


Page 

40 


Steady  State  Problem  Solution  . 

2  Steady  State  Error  at  Finite  Difference  Mesh 


Points .  41 

3  Steady  State  Error  at  99  Points  .  43 

4  Transient  Problem  Exact  Solution  at  Three 

Dimensionless  Times  .  49 

5  Transient  Error  for  Time  =  .05 .  50 

6  Transient  Error  for  Time  =  .10 .  51 

/  Transient  Error  for  Time  =  .15 .  52 


iv 


List  of  Tables 


Table  Page 

I  Error  and  Time  for  Collocation  on  the  Steady 

State  Problem .  46 

II  Error  and  Time  for  Collocation  in  the 

Transient  Problem  .  55 

III  Error  and  Time  for  Galerkin  on  the  Steady 

State  Problem .  68 

IV  Error  and  Time  for  Least  Squares  on  the 

Steady  State  Problem  .  69 

V  Error  and  Time  for  Finite  Difference  Method 

on  the  Steady  State  Problem .  70 

VI  Error  and  Time  for  Galerkin  in  the  Transient 

Problem .  71 

VII  Error  and  Time  for  the  Explicit  Finite 

Difference  Method  in  the  Transient  Problem 

T=.  05 .  72 

VIII  Error  and  Time  for  the  Explicit  Finite 

Difference  Method  in  the  Transient  Problem 

T=.  10 .  72 

IX  Error  and  Time  for  the  Explicit  Finite 

Difference  Method  in  the  Transient  Problem 

T=.  15 .  73 

% 


V 


AFIT/ GNE/PH/8 1-7 


Abstract 

The  method  of  weighted  residuals  applied  on  the  whole 
domain  for  steady  state  and  transient  heat  generation 
problems  was  compared  to  finite  difference  methods.  The 
comparison  consisted  of  the  maximum  absolute  error  from 
the  exact  solution  and  computer  time  required  for  solution. 

One  steady  state  and  one  transient  heat  generation 
problem  were  solved  by  collocation  and  Galerkin  weighted 
residual  methods  and  finite  differences.  The  least  squares 
weighted  residual  method  was  also  used  for  the  steady  state 
problem.  Both  problems  were  one  dimensional  and  had 
Dirichlet  boundary  conditions.  Integrals  for  weighted 
residual  methods  were  evaluated  analytically  to  produce 
recursion  relations.  The  transient  problem  was  solved  by 
the  reduction  to  ordinary  differential  equations  method  for 
weighted  residuals. 

The  Galerkin  method  was  fastest  to  a  given  accuracy 
for  both  problems  evaluated.  The  accuracy  of  Galerkin  and 
other  weighted  residual  methods  was  greater  than  finite 
differences  after  a  point  at  low  solution  accuracy.  This 
crossover  point  was  typically  two  to  three  digits  of  accu¬ 
racy.  The  polynomial  trial  functions  used  for  weighted 
residual  solutions  exhibited  a  numerical  instability  for 
solutions  of  10  terms  and  over  increasing  the  maximum  abso¬ 
lute  error.  Orthogonal  collocation  and  weighted  residuals 
on  finite  elements  were  recommended  as  alternate  methods. 


INVESTIGATION  OF  THE  NUMERICAL  METHODS  OF 
FINITE  DIFFERENCES  AND  WEIGHTED  RESIDUALS 
FOR  SOLUTION  OF  THE  HEAT  EQUATION 


I  Introduction 

The  method  of  finite  differences  has  long  been  used  for 
numerical  solution  of  the  heat  equation  and  other  differen¬ 
tial  equations  similar  to  it.  One  example  is  the  diffusion 
equation  of  nuclear  engineering.  Another  method,  the  method 
of  weighted  residuals  used  with  finite  elements,  has  pro¬ 
duced  results  of  equal  accuracy  at  a  fraction  of  the  compu¬ 
ter  time  (Ref  6:323-350).  These  results  were  obtained  for 
problems  of  structural  stress  which  are  similar  to  the 
steady  state  heat  equation.  This  paper  will  investigate  the 
method  of  weighted  residuals  (MWR)  applied  over  the  whole 
domain  of  the  problem  to  see  if  it  has  similar  advantages 
over  finite  differences. 

The  objective  of  this  study  is  to  compare  and  evaluate 
the  finite  difference  method  and  the  method  of  weighted 
residuals  on  the  whole  domain  for  the  steady  state  and 
transient  heat  equations.  The  basis  of  the  comparison  will 
be  accuracy  of  solution  and  computer  time  required.  For 
this  comparison  the  MWR  will  be  used  on  the  whole  domain  of 
interest  and  not  finite  elements.  One  steady  state  and  one 
transient  problem  are  solved  using  both  methods  and  the 
results  compared.  Both  problems  are  in  one  space  dimension 
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of  a  cartesian  coordinate  system  with  Dirichlet  boundary 
conditions . 

The  paper  consists  of  chapters  on  the  theory  of  MWR  and 
finite  differences,  the  problems  chosen,  results,  and  con¬ 
clusions.  The  theory  chapter  discusses  the  MWR  for  steady 

| 

state  and  transient  problems  and  finite  differences.  In  the 
problems  chapter,  the  problems  are  discussed  and  the  MWR  and 
finite  difference  solutions  are  developed.  The  results  and 
comparisons  chapter  lists  the  results  of  comparison  and 
problems  encountered.  Last,  the  conclusion  chapter  gives 

the  overall  conclusions  of  the  investigation.  i 
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II  Theory 


The  MWR  used  on  the  whole  domain  of  interest  constructs 
a  series  function  which  approximately  satisfies  the  differ¬ 
ential  equation  and  boundary  conditions  over  the  domain. 

The  finite  difference  method,  on  the  other  hand,  produces  a 
table  of  values  which  approximate  the  solution  to  the  dif¬ 
ferential  equation  and  boundary  conditions  at  the  mesh 
points.  Values  of  the  solution  between  mesh  points  must  be 
obtained  by  interpolation  for  finite  differences,  but  are 
obtainable  directly  from  the  MWR  series  solution  for  any 
point  in  the  domain. 

The  comparison  between  methods  will  be  made  by  measuring 
the  time  to  attain  various  accuracies.  Accuracy  is  increased 
by  finer  mesh  spacing  for  finite  differences  and  more  terms 
in  the  series  solution  for  MWR.  This  chapter  will  cover  the 
MWR  and  finite  difference  methods  to  be  compared. 

The  Method  of  Weighted  Residuals 

The  basic  idea  of  the  MWR  is  a  series  of  complete  and 
linearly  independent  trial  functions  which  is  made  to  satisfy 
the  differential  equation  and  boundary  conditions  (Ref  3:35). 
The  series  is  substituted  into  the  equation  and  boundary 
conditions  and  the  error  required  to  vanish  in  an  average 
sense. 

Take  the  case  of  general  boundary  conditions  for  the 
steady  state  problem: 
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LU(x) 


f (x)  in  Region  D 


(1) 


KU(x)  =  ^<s)  on  Boundary  S  (2) 

where  L  is  the  differential  operator  and  K  is  the  boundary 
condition  operator.  The  solution,  U  (x) ,  is  approximated  as 

N 

U  (x)  =  E  a.  <J>.  (x)  (3) 

i=l  1  1 

where  a^  are  undetermined  constants  and  <j;^  (xj  are  the  trial 
functions  which  make  up  the  trial  solution  expansion  (Ref  8 : 
258-261,  3:8).  Substituting  UN(x)  for  U(x)  in  Eqs.  (1)  and 
(2)  yields  the  residual,  R: 


LUn(x) 

-  f (x)  =  R 

(4) 

KUn(x) 

-  i|»(s)  =  R 

(5) 

Since  UN(x)  is  just  an  approximation,  the  equations  may  not 
balance.  The  error  or  residual  of  UN (x)  instead  of  U(x)  in 
Eqs.  (4)  and  (5)  is  integrated  with  a  weight  function  to 
minimize  that  error: 


^wj(LUN(x) 

jH'*1 


f  (x))  dD 
ip  (s))  dS 


j=l, . . . N 

(6) 

j=l , . .  .  N 

(7) 

The  weight  functions,  w^ ,  may  be  determinea  in  several  ways 
and  will  be  discussed  later.  There  must  be  N  weight  func¬ 
tions  to  solve  for  the  N  undetermined  constants. 
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The  undetermined  constants  in  U.,  (x)  are  found  by  solv- 

N 

ing  Eqs .  (6)  and  (7)  simultaneously.  Substituting  the 
expansion  for  U„ (x)  of  Eq.  (3)  into  Eqs.  (6)  and  (7)  yields 


for  j  =  1, . . .N 

N  r  _ 

E  a-.  Jw.4,  (x)  dD 
i=l  1  D  3  1 

N  r 

E  a.  Jw.  K4> .  (x)  dS 

i=l  1  s  3  1 


-  f  w.  f  (x)  dD  =  0  (8) 

D  3 


-  y*w.tMs)dS  =  0  (9) 

S  3 


Equating  Eqs.  (8)  and  (9)  yields  for  j  =  1,...N 


N 

E 

i=l 


ai 


f w.Ltfi.  (x)dD  -  /w.K*.  (x)dS 
D  3  1  S  3  1 

=  f  w^f(x)dD  -  ^/w  •  ^  (s )  dS 
U  J  C  J 


(10) 


Thus  Eq.  (10)  becomes  an  N  by  N  matrix  problem  of  simultane¬ 
ous  equations  in  a^ 

A  a  =  b  (11) 


where  a  is  the  vector  of  arbitrary  constants  a^,  where  b  is 
the  vector  of  integrals  not  involving  a^  for  j  =  1,...N 


b  =  /w,f(x)dD  -  yw.^(s)dS 
D  3  S  3 

and  where  A  is  the  coefficient  matrix  for  j  =  1,...N 


(12) 


A  = 


N 

E 

i=l 


fvr.LQ.  (x)dD  -  /w'Ki).  (x)  dS 
n  J  1  <5  3  1 


(13) 
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Note  that  an  underlined  capital  letter  symbolizes  a  square 
matrix  and  an  underlined  lower  case  letter  symbolizes  a 
column  matrix  or  vector.  Once  the  system  of  simultaneous 
equations  has  been  solved  for  the  a^,  the  approximate  solu¬ 
tion  has  been  determined  for  the  whole  domain. 

Types  of  Method  Application.  There  are  three  ways  in 
which  MWR  is  generally  applied.  They  are  boundary  methods, 
interior  methods,  and  mixed  methods  (Ref  3:11). 

In  boundary  methods  the  trial  functions  satisfy  the 
differential  equation  exactly,  but  not  the  boundary  condi¬ 
tions.  MWR  is  then  applied  on  the  boundary  conditions  only. 
Then  Eq.  (4)  for  differential  equation  residual  is  satisfied 
exactly  and  those  integrals  in  Eqs.  (6)  and  (8)  total  zero. 
Last,  Eq.  (10)  which  specifies  the  matrix  equations  to 
solve  for  unknown  constants  becomes 


(14) 


For  interior  methods  the  trial  functions  satisfy  the 
boundary  conditions  so  MWR  is  applied  on  the  differential 
equation  only.  Then  as  with  boundary  methods  Eqs.  (7)  and 
(9)  are  no  longer  needed  and  Eq.  (10)  becomes  for  j  =  1,...N 


N 

E 

i=l 


i.  /w.L4>.(x)dD  =  /  w.f(x) 

1  |*D  31  V 


The  problems  studied  in  this  paper  are  solved  with  interior 
MWR. 
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Last  in  mixed  methods  the  trial  functions  may  satisfy 
neither  the  differential  equation  nor  boundary  conditions. 
The  trial  functions  may  satisfy  some  parts  of  either.  A 
boundary  condition  that  is  usually  satisfied  is  the  value  of 
the  function  on  the  boundary  (Ref  3:28-30).  Eq.  (10) 
remains  the  same  for  mixed  MWR. 

Types  of  Weighting  Criteria.  Several  types  of  criteria 
for  minimizing  the  residual  exist.  Each  criteria  has  an 
associated  set  of  weight  functions  which  are  used  to  inte¬ 
grate  the  residuals  in  MWR. 

The  weight  function  in  the  collocation  method  is  the 
Dirac  delta  function  (Ref  2:148): 

Wj  =  6  (x  -  x j )  (16) 

The  delta  function  f  Drees  the  value  of  the  residual  integral 
of  Eqs.  (6)  and  (7)  to  equal  zero  at  the  collocation 
point  specified  in  the  delta  function.  For  a  trial 
function  expansion  of  N  terms,  N  points  are  needed  to  gen¬ 
erate  N  equations  for  simultaneous  solution. 

Choice  of  the  collocation  points  can  have  a  definite 
effect  on  the  accuracy  of  the  solution  obtained  for  expan¬ 
sions  with  a  low  number  of  terms.  A  usual  practice  is  to 
choose  them  evenly  spaced  through  the  region  or  boundary. 

At  any  event  as  the  residual  is  made  zero  at  more  and  more 
points,  it  presumably  approaches  zero  throughout  the  region 
or  boundary  (Ref  3:9).  Orthogonal  collocation  uses 
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orthogonal  polynomials  as  trial  functions  and  chooses  the 
collocation  points  as  the  roots  to  the  polynomials  to 
improve  the  accuracy  of  low  order  expansions  (Ref  3:97-98). 
Also  a  method  called  least  squares-collocation  uses  more 
collocation  points  than  constants  to  be  determined.  The 
residuals  are  squared,  added,  and  minimized  with  respect  to 
the  constants  (Ref  3:26-27).  Least  squares-collocation  on 
the  boundary  has  been  applied  as  a  method  for  steady-state 
and  transient  heat  problems  with  arbitrary  geometry  and 
boundary  conditions  in  two  and  three  dimensions  (Ref  16:103- 
108)  . 

A  second  weighting  criteria  called  the  subdomain  method 
integrates  the  residual  over  N  subintervals  of  the  domain. 
Then  the  weight  function  is  (Ref  2:149) 


{1  in  subdomain  i 
0  outside  of  subdomain  i 


(17) 


The  Galerkin  method  corresponds  to  a  third  criteria  for 
minimizing  the  residual.  In  the  Galerkin  method  the  trial 
functions  in  the  approximate  solution  expansion  of  Eq.  (4) 
are  used  as  weight  functions: 

Wj  =  (18) 

This  amounts  to  making  the  weight  functions  orthogonal  to 
the  residual  (Ref  2:149). 

A  fourth  method  is  least  squares.  Here  the  integral  of 
the  square  of  the  residual  is  minimized  with  respect  to  the 


solution  constants  (Ref  2:150) 

8 


JL  / 

3ai  D 


R2  dx  =  2 


/ 


3R 
D  8ai 


R  dx  = 


(19) 


where  R  is  the  residual,  D  represents  the  domain,  and  a^  are 
the  solution  constants  of  Eq.  (4) .  The  weight  function  for 
least  squares  is  then  the  derivative  of  the  residual  with 
respect  to  one  of  the  solution  constants: 


3R 

wi  3a, 


(20) 


This  weight  function  then  becomes  the  a^  term  of  the 
residual . 

Other  weight  functions  than  these  are  used.  In  fact 
the  weight  function  can  be  any  complete  set  or  N  members  of 
it  (Ref  3:11) . 


The  criteria  used  for  the  problems  in  this  paper  will 
be  collocation,  Galerkin,  and  least  squares.  These  were 
the  criteria  used  in  the  paper  comparing  FDM  with  the  MWR 
applied  in  finite  elements  as  mentioned  before  (Ref  6:323- 
350)  . 


Trial  Functions .  The  trial  functions  provide  much  of 
the  power  of  MWR  since  they  incorporate  known  information 
into  the  solution.  Trial  functions  can  incorporate  the 
general  symmetry  of  a  problem  or  satisfy  some  of  the  boundary 
conditions  (Ref  3:35).  If  the  boundary  conditions  in  a 
problem  are  of  the  first  kind,  such  as  UQ(x)  =  f  (x)  =  i|»(s), 
then  a  convenient  form  of  trial  solution  is 


Vx) 


£  (x) 
o 


N 

+  E  a.  .  (x) 
i=l  1 


(21) 
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where  the  trial  functions  <j>^(x)  are  zero  on  the  boundary  and 
fQ(x)  reduces  to  if>(s)  on  the  boundary.  Then  a  problem  like 
Eq.  (1)  and  (2)  can  be  transformed  into  an  interior  MWR 
problem.  Also  fQ(x)  which  represents  the  boundary  condition 
can  be  eliminated  from  the  trial  solution  by  applying  the 
differential  operator  of  Eq.  (1)  to  the  trial  solution  term 
by  term.  Then  Eq.  (1)  becomes 


LU  (x)  =  Lf  (x)  +  LU  (x)  =  f(x)  (22) 

O  ON 


where  L  is  the  differential  operator,  f (x)  is  defined  as  in 
Eq.  (21),  and  UN(x)  is  the  new  trial  solution 


N 

U  (x)  =  U  (x)  -  f  (x)  =  Z  a.<{>.  (x) 

N  o  o  i=l  ^  i 


(23) 


Then  Eqs.  (1)  and  (2)  become  a  problem  with  a  new  non- 
homogeneous  part 


LU(x)  =  f (x)  -  LfQ (x) 
and  a  new  boundary  condition  (Ref  3:30) 

U (x)  =  0 


(24) 


(25) 


The  set  of  functions  chosen  as  trial  functions  ij>^(x) 
as  in  Eq.  (23)  must  be  complete  and  linearly  independent  to 
represent  the  solution  to  a  boundary  value  problem  (Ref  3: 
35).  One  such  set  is  the  polynomials.  Linearly  independent 
and  continuous  polynomials  have  been  proven  complete  for 
homogeneous  and  nonhomogeneous  steady  state  heat  equations 
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in  one  and  two  dimensions  (Ref  3  :  355-356 , 359 , Ref  8:263-265, 
273,277).  In  this  case  the  solutions  evolved  are  uniformly 
convergent  to  the  true  solution.  Then  we  can  choose  an 
e  >  0  and 


I  ^true  ^  -  Vx>l  <  £  (26> 

where  U.  (x)  is  the  true  solution  and  U>T  (x)  is  the  trial 
true  N 

solution  as  in  Eq.  (23)  .  The  property  of  uniform  convergence 
provides  a  test  of  accuracy  since  the  absolute  value  of  the 
error  should  decrease  for  trial  solutions  which  are  expan¬ 
sions  of  more  terms  of  the  set  of  trial  functions.  As  more 
terms  are  taken  the  trial  solution  should  represent  the 
true  solution  more  accurately.  However,  the  exact  choice  of 
the  trial  functions  will  influence  accuracy  in  low  order 
expansions  and  affect  the  rate  of  convergence  (Ref  3: 

34-36)  . 

Trial  functions  for  time  dependent  problems  such  as 
the  transient  heat  equation  can  be  obtained  using  trial 
functions  in  the  space  variables  which  satisfy  the  boundary 
conditions  multiplied  by  unknown  functions  of  time: 

N 

U„(x,t)  =  Z  A.  (t)<f>.(x)  (27) 

N  i=l  1  1 

where  A^(t)  is  the  unknown  function  of  time  (Ref  3:36).  The 
time  functions  are  found  by  applying  MWR  and  the  initial 
conditions  as  illustrated  in  the  next  section  (Ref  3:44-45). 

Transient  Problems .  The  MWR  solution  methods  developed 
earlier  in  this  chapter  were  for  steady  state  problems  only. 
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For  transient  problems  using  the  trial  functions  of  Eq.  (27) 
the  solution  process  is  different.  The  dependence  of  the 
trial  functions  on  time  as  well  as  a  multiplying  constant 
must  be  determined.  Applying  the  MWR  to  the  boundary  condi¬ 
tions  and  differential  equation  determines  the  time  depend¬ 
ence.  Then  applying  the  MWR  to  the  initial  conditions  with 
time  equal  to  zero  in  the  trial  solution  determines  the 
constant  multipliers  (Ref  3:44-45).  Such  a  treatment  is 
called  reduction  to  ordinary  differential  equations. 

The  determination  of  the  time  dependence  of  the  trial 
solution  for  the  transient  heat  equation  requires  solving 
simultaneous  ordinary  differential  equations  in  time.  The 
transient  heat  equation  in  dimensionless  form  with  boundary 
and  initial  conditions  in  addition  is 

V2UU,t)  +  f(x,t)  =  (^-t}  (28) 

where  U(x,t)  is  the  dimensionless  temperature,  f(x,t)  is  a 

2  . 

heat  generation  term,  and  V  is  the  Laplacian.  Note  that 
heat  generation  depending  on  temperature  is  not  considered 
here.  If  the  trial  function  of  Eq.  (27)  which  satisfies 
the  boundary  conditions  is  substituted  into  Eq.  (28) ,  the 
residual,  R,  is 

Nf  2-1  -  N  f  -  5Ai(t)l 

z  [Ai(t)v  4>i(x)J  +  -  .£1[<fri(*)-  at  "J  =  R  (29) 

Applying  a  MWR  weight  criteria  to  this  residual  yields  for 
j  =  1 , . . .  N 
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N  r  2  -  n  r 

2  A.  (t)  J  W.V  <J>  -  (x)dD  =  I  A.  (t)  /  w.<j> .  (x)dD 
i=l  1  3  1  i=l  3t  1  •'d  3  1 


-/ 


w.f  (x,  t)dD 


(30) 


where  w^  is  the  weight  function  for  the  criteria  used.  Eg. 
(30)  can  be  expressed  in  matrix  form  as 


G  cl 


b  (t) 


(31) 


where  G  and  H  are  the  coefficient  matrices  representing  the 
integrals  in  Eq.  (30),  where  a  is  the  vector  of  A^(t),  and 
where  b(t)  is  the  vector  representing  the  last  integral  of 
Eq.  (30)  (Ref  4:735-736).  If  the  inverse  of  H  can  be  found, 
Eq.  (31)  can  be  reduced  to 

^  a  =  H_1  G  a  +  H_1  b ( t)  (32) 

where  H  ^  is  the  inverse  of  H.  Eq.  (32)  is  a  system  of 
simultaneous,  linear,  ordinary  differential  equations  for 
i  =  1  ,N 

dA. (t) 

~5t —  =  ailAl(t)  +  •*•  +  aiNAN(t)  +  fi(t)  (33) 


where  A^ (t)  is  an  element  of  vector  a ,  where  a|^  . . .  a|N  are 
elements  of  matrix  H*"^  G,  and  f ^  (t)  are  elements  of  vector 
H_1  b (t)  (Ref  9:218-219). 

The  f^(t)  term  in  Eq.  (33)  makes  it  a  nonhomogeneous 
system  of  first  order  differential  equations.  If  the 
heat  generation  term  f(x,t)  in  Eq,  (28)  is  a  function 
of  position  only,  the  problem  may  be  separated  into  a 
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steady  state  and  a  homogeneous  transient  problem.  Substitut¬ 
ing  a  new  variable  for  temperature  in  Eq.  (28)  and  its  bound¬ 
ary  and  initial  conditions  can  separate  the  problem: 

T(x,t)  =  U(x,t)  +  V (x)  (34) 

where  T(x,t)  is  the  new  temperature  variable,  U(x,t)  is  the 
temperature' for  the  transient  problem,  and  V(x)  is  the  tem¬ 
perature  for  the  steady  state  problem  (Ref  11:152-153).  The 
steady-state  problem  can  be  solved  by  MWR  developed  earlier 
in  this  chapter.  The  homogeneous  transient  problem  is  then 
solved  by  the  MWR  reduction  to  ordinary  differential  equa¬ 
tions  outlined  in  this  section. 

For  the  homogeneous  transient  problem,  Eq.  (28)  reduces 
to 


V2U (x, t) 

au(x,t) 

3t 

(35) 

U  (s ,  t) 

=  4>  (s) 

(36) 

U(x,0) 

=  g  (X) 

(37) 

where  Eq.  (36)  represents  boundary  conditions  which  do  not 
vary  with  time  and  Eq.  (37)  is  the  initial  condition.  The 
transient  problem  considered  in  this  paper  is  described  by 
Eq.  (35) ,  (36)  and  (37) .  Again  .substituting  the  trial 
function  of  Eq.  (27)  which  satisfies  the  boundary  conditions 
into  Eq.  (35)  and  applying  MWR  leads  to  a  statement  similar 
to  Eq.  (30)  for  j  =  1,...N: 
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{ 


N 

Z 

i=l 


Ai(t) 


i(x)dD 


N 

Z 

i=l 


Jt  Ai<t> 


^  (x)  dD 


(38) 


where  w^  is  the  weight  function  and  A^  (t)  the  unspecified 
function  in  time.  Equation  (38)  expressed  in  matrix  form  is 

da 

—  —  =  —  dt  (39) 

where  G  and  H  are  the  coefficient  matrices  as  before  and  a 
is  the  vector  of  A^(t).  Multiplying  by  the  inverse  of  H 
again  yields  a  system  of  ordinary  differential  equations  in 
time  which  are  homogeneous  in  this  case: 

dA  (t) 

-it—  -  ahAi(t>  +  •••  +  •ihV41  1401 


where  a!,  ...  a'..,  are  the  elements  of  the  coefficient  matrix 
xl  xN 

H  3  G  as  before. 

The  general  solution  to  the  problem  of  Eq.  (40)  is 
(Ref  9:220-222,294) 

Alt  A^t  A^t 

Ai(t)  =  Claile  +  C2ai2e  +  •••  +  cNaiNe  (41) 


where  A^  ...  AN  are  the  eigenvalues  of  the  system  of  equa¬ 
tions  given  in  Eq.  (39)  or  more  specifically  the  matrix 
H  1  G.  Also  the  a.,  .  ...  a„ .  are  the  components  of  the 

-  -  id  n;i 

linearly  independent  eigenvectors  of  matrix  H  G  and 

...  CN  are  the  unspecified  multipliers  which  multiply  all 
elements  of  that  eigenvector.  Eq.  (39)  is  then  an  eigen¬ 


value  problem  since  the  derivative  of  an  eigenvector  a. .e 
At  13 

j  Q  .  OQA \  . 


is  a . . A  .e 
J-D  D 


(Ref  9:294)  : 
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G  a  =  A  H  a 


(42) 


Forming  the  solution  for  (t)  is  then  a  two-step  pro¬ 
cess.  First  the  eigenvalues  and  eigenvectors  of  the  system 
of  ordinary  differential  equations  describing  the  problem 
are  found  to  determine  the  time  dependence  of  exponential 
decay  of  the  transient  solution.  The  last  step  of 
determining  the  constant  multipliers,  ...  C N,  is 

analogous  to  the  MWR  for  steady  state  problems.  The  trial 
solution  for  time  zero  is  substituted  into  the  initial  con¬ 
dition,  Eq.  (37) ,  to  form  a  residual.  The  residual  is  then 
minimized  with- a  weight  function  to  yield  simultaneous 
equations  for  ...  CN  with  j  =  1,...N: 


N  f 

l  A. (0) J w  .  (j>  . (x)dD 
i=l  1  3 


/ 

•'d 


w^g (x)dD 


(43) 


where  A^(0)  is  the  time  function  of  Eq.  (27)  evaluated  at 
the  initial  time,  where  g(x)  is  the  initial  temperature  dis¬ 
tribution,  Wj  the  weight  function,  and  (x)  the  trial  func¬ 
tions.  Since  A^(t)  is  a  sum  of  exponentials  in  time,  for 
the  initial  time  zero  they  take  a  specific  value  of  one. 

Each  A^(0)  then  becomes  a  sum  of  the  eigenvector  multipliers, 
...  CN  and  the  summation  can  be  changed  to  add  the 
terms  from  all  A^(0)  and  the  same  for  the  other  multipliers. 
The  matrix  system  then  reduces  to 

J  c  =  d  (44) 
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where  J  is  the  coefficient  matrix  of  multipliers,  c  is  the 

multiplier  vector,  and  d  is  the  vector  of  initial  values. 

Once  the  simultaneous  equations  are  solved,  the  (t)  and 

MWR  solution  is  finished  (Ref  9:294). 

Only  two  of  the  three  MWR  criteria  used  for  steady 

state  problems  are  applicable  to  this  transient  method. 

Collocation  and  Galerkin  both  result  in  the  eigenvalue 

problem  of  Eq.  (42) .  However,  least  square  does  not  since 

2 

the  resulting  equations  have  factors  of  X  as  well  as  X. 

Only  the  collocation  and  Galerkin  criterias  will  be  used 
for  the  transient  problem  (Ref  2:313-314). 

Numerical  Solution  for  MWR 

This  section  will  outline  the  solution  algorithms  for 
the  steady  state  and  transient  problems. 

Steady  State.  The  numerical  solution  of  MWR  on  the 
whole  domain  of  interest  has  been  largely  limited  to  collo¬ 
cation  methods.  These  methods  include  least  squares-collo- 
cation  and  orthogonal  collocation  mentioned  earlier  under 
the  section  on  types  of  criteria.  Collocation  offers  an 
advantage  since  no  integrals  need  to  be  evaluated.  Only 
the  values  of  residuals  at  collocation  points  need  be  eval¬ 
uated. 

Collocation,  Galerkin,  and  least  square.,  criterias 
were  used  to  compare  MWR  on  finite  elements  with  finite 
differences  in  an  article  by  Houstis  et.  al.  (Ref  t>:323-350). 
In  order  to  compare  the  results  in  this  paper  with  the 
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results  of  Houstis,  all  three  criteria  will  be  used  for  the 
steady  state  problem. 

Evaluation  of  the  integrals  for  the  Galerkin  and 
least  squares  criteria  presents  another  difficulty.  A  form 
of  numerical  integration  could  be  used  to  evaluate  the 
integrals  over  the  domain  of  interest.  However,  since  the 
trial  functions  used  are  simple  polynomials  and  the  integra¬ 
tion  on  one  dimension  only  for  the  problems  in  this  paper, 
the  integrals  can  be  easily  evaluated  analytically  to 
develop  a  recursion  relation  for  the  matrix  elements.  The 
recursion  relation  then  needs  to  be  evaluated  only  once  for 
each  matrix  element  where  the  integrand  would  need  to  be 
evaluated  several  times  for  numerical  integration. 

Another  choice  to  be  made  is  the  choice  of  a  matrix 
equation  solver.  Since  every  element  of  the  MWR  coefficient 
matrix  will  have  a  non-zero  value,  the  best  choice  appears 
to  be  Gaussian  elimination  as  a  matrix  equation  solver. 
Iterative  methods  can  be  faster  than  Gaussian  elimination 
for  sparse  matrices,  but  the  rate  of  convergence  is  uncer¬ 
tain.  The  Gaussian  elimination  will  be  used  for  both  MWR 
and  finite  difference  solutions. 

Then  the  algorithm  for  solution  of  the  steady  state 
problem  for  MWR  consists  of  two  .steps.  First,  the  solution 
matrices  are  formed  from  recursion  relations.  Then  the 
matrix  systems  are  solved  to  obtain  the  defining  constants 
for  MWR. 
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After  the  solution  is  obtained  it  is  compared  against 
the  true  analytical  solution  to  determine  its  accuracy. 

Two  measures  of  accuracy  will  be  used.  The  first  is  abso¬ 
lute  accuracy  which  should  decrease  as  the  MWR  solution  con¬ 
tains  more  terms.  The  second  is  fractional  or  percentage 
accuracy  which  is  most  of  interest  in  a  practical  sense. 

The  fractional  error  determines  the  number  of  significant 
figures  in  the  answer. 

After  the  solutions  are  obtained  and  the  accuracy  is 
checked  for  solutions  with  different  numbers  of  terms,  the 
time  to  compute  the  solution  is  measured.  The  central  pro¬ 
cessor  time  of  the  Control  Data  Corporation  Cyber  74/Cyber 
750  computer  system  with  NOS/BE  operating  system  was  used 
to  measure  the  solution  formation  time.  Since  the  central 
processor  time  output  is  only  accurate  to  .01  seconds,  the 
solution  for  a  given  number  of  terms  must  be  repeated  several 
times  for  accurate  measurement  (Ref  5:8-9). 

Transient.  Numerical  solution  of  the  transient  MWR 
problem  involves  solution  of  an  eigenvalue  problem  and  a 
system  of  simultaneous  linear  equations.  Solution  of  the 
simultaneous  equations  will  be  by  Gaussian  elimination  as 
for  the  steady  state  problem.  The  eigenvalue  problem  will 
be  solved  in  its  general  form. 

The  eigenvalue  problem  is  shown  in  Eq.  (42)  which  is 
repeated  as  Ea.  (45)  here 

G  a  =  X  H  a  (45) 
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where  H  is  the  coefficient  matrix  for  integrals  which  mini¬ 
mize  the  residual  terms  involving  the  time  derivative,  where 
G  is  the  coefficient  matrix  for  the  integrals  which  minimize 
the  space  derivative,  where  a  is  the  vector  of  unknown  time 
functions,  A^(t),  in  the  trial  solution,  and  where  X  is  the 
eigenvalue.  The  vector  a  is  the  set  of  the  N  linearly 
independent  eigenvectors  of  the  system  and  A  represents  the 
N  eigenvalues  associated  with  the  eigenvectors.  The  total 
solution  for  A^ (t)  is  shown  in  Eq.  (41).  A  subroutine  which 
solves  the  eigenvalue  problem  without  inverting  the 
H  matrix  was  selected  from  the  IMSL  program  library.  This 
subroutine,  EIGZF,  allows  the  G  and  H  matrices  to  be  input 
directly  without  inversion  (Ref  10 : 241-256 , 7 :EIGZFl-5) . 

The  algorithm  for  computing  the  transient  MWR  solution 
involves  four  steps.  First  the  G  and  H  matrices  are  computed 
from  Eq.  (38) .  Next  the  eigenvalues  and  eigenvectors  are 
found  by  EIGZF  subroutine.  Then  the  coefficient  matrices  for 
calculation  of  the  constant  multipliers  are  evaluated  as 
indicated  in  Eq.  (43)  and  (44).  Last  the  system  is  solved 
for  the  multipliers. 

As  for  the  steady  state  problem,  the  accuracy  of  solu¬ 
tion  and  formation  time  of  solution  must  be  measured.  The 
one  MWRvSolution  is  good  for  all  times  due  to  the  exponen¬ 
tial  time  dependence  of  the  time  function  A^ (t)  of  Eq.  (27). 
Time  for  solution  must  be  measured  for  different  numbers 
of  terms  in  the  MWR  solution.  Then  accuracy  at  different 
times  is  measured  against  the  analytical  solution. 
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The  Finite  Difference  Method 


The  finite  difference  method  approximates  the  deriva¬ 
tives  of  a  differential  equation  by  the  ratios  of  differ¬ 
ences  between  points  in  the  domain  of  the  problem.  The 
system  of  equations  that  results  can  be  solved  for  the  values 
of  the  function  in  the  differential  equation  at  the  mesh 
points.  Since  the  finite  differences  only  approximate  the 
derivatives,  the  values  for  the  function  have  some  error 
based  on  the  mesh  spacing. 

The  two  derivatives  in  the  steady  state  and  transient 
heat  equations  will  be  approximated  by  central  and  forward 
differences.  The  two  derivatives  are  the  second  derivative 
with  respect  to  the  space  variable  and  the  first  derivative 
with  respect  to  time: 

32U (x, t)  =  62U  (x , t)  =  U (x+h , t )  -  2U (x , t)  +  U(x-h,t) 

(46) 


3x 


3U  (x ,  t)  _  AU (x , t)  _  U  (x,  t+k)  -  U  (x ,  t) 
9t  “  k  k 


(47) 


where  Eq.  (46)  shows  the  second  derivative  with  respect  to 
x,  Eq.  (47)  shows  the  first  derivative  with  respect  to  time, 
where  U(x,t)  is  the  function  differentiated,  where  h  and  k 
are  the  distances  between  mesh  points  in  x  and  time  respec¬ 
tively,  where  6  denotes  a  central  difference  taken  about  the 
point  (x,t)  where  the  derivative  is  approximated,  and  where 
A  denotes  a  forward  difference  taken  from  the  point  (x,t) 
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(Ref  15:59-60).  Eq.  (46)  and  (47)  are  also  valid  for  ordi¬ 
nary  derivatives. 

For  the  steady  state  problem  the  unknown  value  of  the 
function  at  all  mesh  points  will  be  found  simultaneously. 

The  differential  equation  using  finite  difference  approxi¬ 
mations  is  formulated  for  each  mesh  point.  Then  the  system 
of  equations  is  solved  simultaneously  by  the  same  Gaussian 
elimination  method  used  for  MWR. 

Explicit  methods  will  be  used  for  the  transient  prob¬ 
lem.  In  explicit  methods,  the  differential  equation  with 
finite  difference  approximations  is  used  to  compute  the 
value  at  a  next  point  from  values  at  known  points . 

Since  the  difference  expressions  of  Eqs.  (46)  and  (47) 
are  derived  by  combining  Taylor  series  expansions  of  the 
function  to  be  differentiated,  the  error  in  the  expressions 
can  be  found  by  examining  the  expansion  terms  not  used  in 

the  expression.  The  central  difference  then  has  an  error 

,  2  2  2 
proportional  to  h  or  of  order  h  ,  0(h  ),  where  h  is  the 

mesh  spacing.  The  first  forward  difference  has  an  error 
0(k)  where  k  is  the  mesh  spacing.  The  mesh  spacing  deter¬ 
mines  the  error  then  and  closer  mesh  spacing  will  give  more 
accurate  results.  Also  note  that  for  an  equation  that  uses 
two  difference  expressions,  such  as  the  transient  heat  equa¬ 
tion,  the  error  is  the  sum  of  the  two  errors  (Ref  15:59-60, 
108)  . 
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Ill  The  Problem  Set 


One  steady  state  and  one  transient  heat  problem  are 
solved  by  the  MWR  on  the  whole  domain  and  by  the  finite 
difference  method.  Both  problems  are  in  one  space  dimension 
with  Dirichlet  boundary  conditions.  This  chapter  discusses 
the  problems  and  their  solutions .  To  prevent  any  advantage 
between  methods  all  MWR  and  finite  difference  solutions  were 
programmed  in  Fortran  4  Extended.  The  programs  were  also 
compiled  with  the  same  compilation  option,  option  1. 


Steady  State  Problem 

For  the  steady  state  problem  consider  for  0  <  x  <  1 


d2U (x) 
dx2 


+  U  (x)  +  x  =  0 


(48) 


U  (0)  =  U(l)  =  0  (49) 

where  Eq.  (48)  represents  the  heat  equation,  where  Eq.  (49) 
is  the  homogeneous  Dirichlet  boundary  condition,  where  U(x) 
is  temperature,  and  where  U(x)+x  represents  a  heat  genera¬ 
tion  term.  Removal  of  the  thermal  conductivity  constant 
factor  which  is  normally  shown  multiplying  the  second  deriv¬ 
ative  of  temperature  and  reduction  of  the  equation  to  this 
form  requires  a  very  special  heat  generation  term.  However, 
study  of  the  equation  in  this  form  provides  a  convenient 
closed  form  solution  (Ref  3:269): 


U(x) 


sin (x) 
sin  (1) 


(50) 
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For  the  MWR  solution  to  the  steady  state  problem  a  poly¬ 
nomial  trial  function  which  satisfies  the  boundary  conditions 
and  has  been  proven  complete  for  similar  problems  is  (Ref  3: 
356) 

=  x^(l-x)  =  x^  —  xi+1  (51) 

Then  the  trial  solution  U^fx)  will  be 

N 

U  <x)  =  E  a.  (x1  -  x1  X)  (52) 

N  i=l  1 


where  a^  are  undetermined  constants.  Substituting  Eq.  (52) 
into  Eq.  (48)  yields  the  residual,  R,  of  the  differential 
equation  for  this  interior  method  application  of  MWR: 

R  =  E  a . [i  (i-1 ) xi_2  -  i(i+l)xi"1  +  x1  -  x1+1J  +  x 
i=l  1 

(53 

Simultaneous  equations  can  now  be  constructed  from  the 
residual  by  integration  with  the  weight  functions  of  the 
three  MWR  methods  to  be  evaluated,  for  j  =  1,...N: 


N  f1  i 

E  a .  I  w . [ i (i-1) x 

i=l  3 


-2  .  , • . i \  i-1  ,  i  i+1,  , 

-  i(i+l)x  +  x  -  x  )dx 


■  -jb 


x  dx 


(54) 


where  Wj  i  *•  the  weight  function  and  Eq.  (54)  represents  one 
of  the  N  simultaneous  equations.  Using  the  weight  functions 
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of  Eq.  (15),  (16),  and  (17)  the  integrals  are  evaluated  to 

give  recursion  relations.  The  recursion  relations  evaluated 
for  values  of  i  and  j  provide  the  matrix  elements  for  the 
matrix  problem  A  a  =  b  where  A  is  the  coefficient  matrix 
defined  by  the  left  side  of  Eq.  (54),  b  is  the  vector 
defined  by  the  right  side  of  Eq.  (54) ,  and  where  a  is  the 
vector  of  undetermined  constants,  a^.  The  specific  weight 
function,  w ^ ,  used  defines  the  rows  of  A  and  the  ith  term 
of  the  trial  solution  defines  the  columns.  Once  the  matrix 
problem  has  been  solved  for  a  the  MWR  solution  is  finished. 

Recursion  relations  for  the  collocation  method  are  the 
residuals : 

A(j,i)  =  ili+Dx^-1  -  i(i-l)x^-2  +  x^+1  -  xj  (55) 

b ( j )  =  Xj  (56) 

where  x^  is  the  jth  of  the  N  equally  spaced  collocation 
points,  A ( j , i)  is  an  element  of  the  matrix  A,  and  b(j)  is 
an  element  of  vector  b.  Also  the  recursion  relations  for 
the  Galerkin  method  using  the  trial  functions  as  weight 
functions  are 


A(j  ,i) 


[(i+l)i  +  i(i-D)  .  2  1 

(j+i)  (j+i+2)  "  ( j+i+3) 


i(i-l)  i  (i+1)  +  1 
( j+i-1)  ”  ( j+i+l) 


(57) 


b  ( j ) 


(j+2Mj+3) 


(58) 
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Last,  the  recursion  relations  for  the  least  squares  method 
with  terms  of  the  residual  as  weight  functions  are 


A(j  ,i) 


j ( j-1) i (i-1)  [j (j-1) i (i+1)  +  i  (i-1) j ( j+1) ] 
(j+i-3)  (i+j-2) 


.  [j  (j-1)  +  j ( j+1) i (i+1)  +  i  (i-1) ] 
(1+3-1) 


[j  ( j-1)  +  j(j+l)  +  i (i+1)  +  i  (i-1) ] 
(i+j) 


b  ( j ) 


,  [i (i+1)  +  1  +  j (j+1) ]  21 

(i+j+1)  (i+j+2)  (i+j+3) 

(59) 


(3+2)  (3+3) 


(60) 


The  A  matrices  resulting  from  the  three  methods  have  all 
non-zero  elements. 

The  finite  difference  solution  is  constructed  using  a 
central  difference  approximation  for  the  second  derivative 
of  temperature  with  respect  to  x.  The  central  difference 
approximation  of  Eq.  (46)  for  temperature  as  a  function  of  x 
only  is 


IZfci 


(61) 


where  the  subscript  j  identifies  the  particular  mesh  point 
of  N  equally  spaced  mesh  points,  where  is  the  second 

derivative  with  respect  to  x,  and  where  h  is  the  spacing 
between  mesh  points.  Then  Eq.  (48),  the  differential  equa¬ 
tion,  becomes  the  difference  equation 
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(62) 


x  . 
J 


where  Eq.  (61)  has  been  substituted  into  Eq.  (48)  and  multi- 
2 

plied  by  -h  .  With  N  mesh  points  the  interval  0  to  1  of  the 
problem  is  divided  into  N+l  spaces  so  h  =  1/(N+1).  Then 
the  matrix  problem  A  a  =  b  for  simultaneous  solution  of  Eq. 
(62)  at  the  mesh  points  can  be  defined  as  for  MWR  where  the 
jth  element  of  a  is  the  value  of  temperature,  LL  ,  at  the  jth 
mesh  point.  The  matrix  A  is  tridiagonal  with  main  diagonal 
elements  (2-h  ) ,  with  elements  of  the  next  diagonals  above 
and  below  the  main  diagonal  of  -1,  and  with  the  other  ele¬ 
ments  zero.  Since  the  Dirichlet  boundary  conditions  are 

2 

homogeneous,  the  elements  of  b  are  h  x^  .  Once  the  matrix 
problem  A  a  =  b  has  been  solved  for  the  value  of  temperature 
at  the  mesh  points,  a,  the  finite  difference  solution  is 
complete. 

The  systems  of  simultaneous  equations  A  a  =  b  for  MWR 
and  finite  differences  will  be  solved  by  Gaussian  elimina¬ 
tion.  A  subroutine,  LEQTlF ,  was  selected  from  the  IMSL 
program  library  to  perform  the  Gaussian  elimination.  LEQTlF 
performs  Gaussian  elimination  with  partial  pivoting,  equi¬ 
libration,  and  the  Crout  algorithm  (Ref  7 : LEQTlF  1-4). 

LEQTlF  also  indi<  ates  when  solution  is  not  possible  due  to 
a  singular  matrix  A  and  tests  the  solution  a  to  insure  it 
agrees  with  matrix  A  to  a  specified  number  of  digits. 


I 
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Transient  Problem 


The  transient  problem  or  initial  value  problem  is  the 
problem  of  a  wall  warming  up.  The  problem  stated  in  dimen¬ 
sionless  form  is  for  0  <  x  <  1 


S2 

-2-5-  U(x,t) 
3x^ 


3 

3 1 


U(x,t) 


(63) 


u(0,t)  =  1  U(l,t)  =  0 


(64) 


U (x, 0)  =  0 


(65) 


where  Eq.  (63)  is  the  partial  differential  equation,  where 
Eq.  (64)  is  the  Dirichlet  boundary  conditions,  where  Eq.  (65) 
is  the  initial  condition,  and  where  U(x,t)  is  the  dimension¬ 
less  temperature.  The  exact  solution  to  the  problem  is 


U  (x,  t)  =  1  -  x 


-  -  l 
77  n-1 


2  2 

sin(mrx)  -n  r  t 

-  e 

n 


(66) 


where  1-x  represents  the  steady  state  solution  that  the 
whole  solution  decays  to  as  time  increases  and  where  the 
balance  of  Eq.  (66)  represents  the  transient  solution 
(Ref  1:93-96) . 

The  MWR  solution  will  be  of  the  reduction  to  ordinary 
differential  equations  form.  The  trial  functions  will  be 
defined  as  in  Eq.  (27)  where  the  space  part  of  the  trial 
function  satisfies  the  boundary  conditions  and  the  undeter¬ 
mined  function  of  time  is  found  to  satisfy  the  differential 
equation  and  initial  conditions.  The  nonhomogeneous 
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boundary  conditions  can  be  satisfied  as  in  Eq.  (21)  by  a 
series  of  trial  functions  which  is  zero  on  both  boundaries 
and  a  leading  function  which  reduces  to  the  value  of  the 
solution  on  both  boundaries.  Then  we  can  choose  a  trial 
function  as 


UN.(x,t) 


N 

1  -  x  +  E  A .  ( t )  (x1 
i=l  1 


i+1 . 
x  ) 


(67) 


where  UM(x,t)  is  the  trial  solution,  where  1-x  satisfies  the 
nonhomogeneous  boundary  conditions,  where  the  spacial  trial 
functions,  (x1  -  x1  ^) ,  satisfy  the  homogeneous  boundary 
conditions,  and  where  A^(t)  is  the  undetermined  function  of 
time.  The  spacial  trial  function,  (x1  -  x1+1) ,  is  the  same 
one  used  for  the  steady  state  problem. 

The  first  step  in  the  transient  MWR  solution  is  forming 
the  matrices  for  the  matrix  eigenvalue  problem,  G  a  =  X  H  a, 
where  G  is  the  coefficient  matrix  that  corresponds  to  the 
left  side  of  Eq.  (63),  H  is  the  matrix  for  the  right  side  of 
Eq.  (63),  a  is  the  vector  of  A^(t)  time  functions,  and  where 
X  is  the  eigenvalue.  First,  form  the  residual  of  the  dif¬ 
ferential  equation  by  substituting  the  trial  function  Eq.  (67) 
into  Eq.  (63) : 


N  .  _  . 

0  +  E  A.  (t)  [i  (i-1)  x1-^  -  (i  +  Dix1  l] 
i=l  1 


N  •  .+1  3A  (t) 

=  0  +  E  (x1  -  x1  +  i] 

i  =  l  3t 


(68) 
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Note  that  the  residual  of  the  leading  term  of  the  trial 
function  is  zero.  Thus,  the  problem  becomes  an 
eigenvalue  problem  as  in  Eq.  (39).  The  form  of  the  time 
function,  A^(t),  for  a  problem  of  this  type  is 

N  X  t 

Ai(tl  -  ^  c„  v  e  (69) 

where  a^k  represents  one  element  of  the  kth  eigenvector  of 
a,  where  represents  the  multiplier  of  the  kth  eigenvector, 
and  X^  is  the  eigenvalue  associated  with  the  kth  eigenvector. 
Since  the  leading  term  of  the  trial  solution  is  the  steady 
state  solution,  change  the  form  of  A^(t)  to  indicate  it  is 
decaying  in  time: 

N  -X  t 

Ai(t)  =  Ck  aik  e  <70> 

Then  the  residual  Eq.  (67)  becomes 

N 

l  (i(i-l)x1-2  -  (i+1) ix1_1]  a 
i=l  K 

N  .  .  , 

=  -  X  Z  [x1  -  x1+1]  a.  (71) 

*  i=l  ~K 


where  a^  is  the  kth  eigenvector  of  a.  Minimizing  the 
residual  with  the  MWR  weight  function  provides  the  matrices 


to  solve  for  X^  and  a^: 


N  f  i- 

I  /  w . [i ( i— 1 ) x1 

i=l  Jq 

N  ,l 


2  -  (i+D  ix^  dx  a^ 


.  [xi'1’1  -  x'*']  dx  a, 

D  & 


30 


where  w ^  is  the  weight  function.  The  left  side  of  Eq .  (72) 
is  matrix  G  and  the  right  side  H  for  the  eigenvalue  problem 
G  a  =  X  Ha# 


The  two  minimizing  weight  criteria  for  the  transient 
problem  are  collocation  and  Galerkin.  The  G  and  H  matrices 
for  collocation  are 


G(j  ,i) 


-2  , 

[(i-l)ix^-2 


i=l 

(i+1) iXj-1] ,  i-2,...N 


(73) 


H  ( j  /  i ) 


(74) 


where  G(j,i)  and  H(j,i)  are  elements  of  matrices  G  and  H  and 
where  x^  is  the  jth  of  N  equally  spaced  collocation  points  on 
the  interval  0  <  x  <  1.  The  Galerkin  matrices  are 


G(j  ,i) 


i(i-l)  _  2i2  (i+1) i 

(j+i-1)  (j+i)  ( j+i+1) 


(75) 


H ( j , i) 


2  11 
(i+j+2)  "  (i+j+1)  (i+j+3) 


(76) 


The  Galerkin  matrices  are  in  general  symmetric,  but  the 

collocation  matrices  are  not.  Least  squares  was  not  used 
2 

since  it  has  X  terms. 

After  the  matrices  are  formed  the  eigenvalue  problem 
is  solved  by  IMSL  library  subroutine  EIGZF  which  solves  the 
real  eigenvalue  problem  G  a  =  X  H  a  for  the  real  or 
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complex  eigenvalues  and  eigenvectors  (Ref  7:EIGZF  1-5).  The 
routine  reduces  G  to  upper  Hessenberg  form  and  B  to  upper 
triangular  form.  Then  the  routine  transforms  G  to  quasi¬ 
upper  triangular  form  (upper  Hessenberg  with  no  two  con¬ 
secutive  subdiagonal  elements  being  nonzero)  while  keeping 
H  in  upper  triangular  form.  Last,  the  routine  calculates 
eigenvalues  and  eigenvectors  through  an  iterative  operation 
(Ref  10:241-255) .  The  IMSL  routine  includes  a  performance 
index  to  evaluate  how  well  the  problem  was  solved  and  indi¬ 
cates  if  the  routine  cannot  converge  to  one  of  the  eigen¬ 
values  (Ref  7:EIGZF  2-3). 

The  possibility  of  complex  eigenvalues  and  eigenvectors 
requires  a  strategy  for  reducing  all  eigenvalues  to  real 
and  the  eigenvectors  to  real  and  linearly  independent.  The 
eigenvalues  of  real  symmetric  matrices  are  real  so  the 
Galerkin  method  should  have  all  real  eigenvalues  and  eigen¬ 
vectors  (Ref  15:24)  .  The  collocation  method  may  have  some 
complex  eigenvalues  and  eigenvectors  due  to  its  nonsymmetric 
matrices.  A  useful  approximation  for  complex  eigenvalues 
is  the  real  part  (Ref  4:739-740).  Since  complex  eigenvalues 
come  in  conjugate  pairs  with  eigenvectors  that  are  also 
conjugates,  taking  the  real  part  of  the  eigenvalue  results 
in  a  repeated  real  eigenvalue.  Two  linearly  independent 
eigenvectors  are  obtained  by  taking  the  real  part  of  the 
eigenvector  for  one  and  the  imaginary  part  for  the  other 
(Ref  9:228-230) . 
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After  the  eigenvalues  and  eigenvectors  have  been  pro¬ 
duced  and  made  real  and  linearly  independent,  the  initial 
conditions  are  used  to  find  the  appropriate  constant  multi¬ 


pliers  for  the  eigenvectors.  The  multipliers  are  the  C 
of  the  Eq.  (70)  form  of  the  time  function,  (t)  .  The 
procedure  in  the  last  two  paragraphs  has  found  the  eigen¬ 
values,  A^,  and  the  eigenvectors,  a_^.  Now  the  initial 
condition  of  Eq.  (65)  is  minimized  with  the  same  weight 
functions  to  provide  simultaneous  equations  for  C^.  Sub¬ 
stituting  Eqs .  (70)  and  (67)  into  Eq.  (65)  and  integrating 

with  a  weight  function,  w^ ,  yields 


/ 


Wj (1-x) dx  + 


•'ft 


N 

Z 


N 


„  ,  i  i+l._ 

Z  (x  -  x  )C.  a  e 

i=l  k=l  1K 


-Ak(0) 


dx 


J 

'  I 


Wj (0)dx 


(77) 


where  N  different  w^  determine  N  different  equations,  where 
the  left  side  is  the  minimization  of  the  trial  solution  of 
Eq.  (67),  and  where  the  right  side  is  the  minimization  of 
the  initial  value.  Note  time  takes  its  initial  value  zero 
in  Eq.  (77) .  The  residual  equation,  Eq.  (77)  ,  can  be 
rearranged  to  provide  simultaneous  equations  for  by 
exchanging  the  order  of  integration  and  summation: 


N 

Z 

k=l 


N 

Jr  Z 
1=1 


(78) 
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Then  the  coefficients  for  are  obtained  by  summing  the 
elements  of  the  kth  eigenvector  multiplied  by  the  integral 


on  the  left  side  of  Eg.  (78)  .  Eq.  (78)  then  becomes  a 
matrix  problem,  M  c  =  m,  where  M  is  the  coefficient  matrix 
of  the  left  side  of  Eq.  (78) ,  where  c  is  the  vector  of  , 
and  where  m  is  the  vector  of  Eq.  (78)  right  sides.  M  c  =  m 
is  solved  by  Gaussian  elimination  for  the  multipliers,  C  , 

K 

and  the  solution  is  then  complete. 

The  integrals  of  Eq.  (78)  are  evaluated  analytically 
to  provide  recursion  relations  for  M  and  m.  For  collocation 
the  recursion  relations  are 


M(j,i)  =  xj  -  Xj+1 

(79) 

m(j)  =  xj  ~  1 

(80) 

where  M(j,i)  is  the  integral  of  the  (j,i)  element  of  M, 
m(j)  is  the  jth  element  of  m,  and  where  x^  is  the  jth  of  N 
equally  spaced  collocation  points.  For  Galerkin  the 
recursion  relations  are 

=  (j+i+1)  "  (j+i+2)  +  ( j+i+3)  (8 


m(j) 


2  11 
(j  +  2)  "  (j+1)  _  (j+3) 


(82) 


The  finite  difference  solution  for  the  transient  prob¬ 
lem  will  be  constructed  using  an  explicit  four  point  method. 
The  differential  equation,  Eq.  (64),  will  be  approximated 
with  a  first  forward  difference  in  time  and  a  central  dif¬ 
ference  in  space: 
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(83) 


where  j  is  the  index  for  mesh  point  location  on  x  with  uni¬ 
form  spacing  of  h^  between  mesh  points  and  where  k  is  the 
index  for  time  mesh  spacing  of  h^  between  time  levels 
(Ref  15:107-108).  Solving  for  unknown  temperature,  U.  ,  .  , 

3  r  K+-L 

gives  the  explicit  equation 


j  ,k+l 


R'Vi.k 


Vi,*’ 


(1  -  2R)U. 


2 

where  R  =  hfc/hx.  Eq.  (84)  is  used  to  solve  for  the  tempera¬ 
ture  at  a  new  time  level  repeatedly  until  the  desired  time 


is  reached.  For  the  explicit  method  to  be  stable,  the  time 

2 

step  size  must  obey  the  stability  criteria  hfc  <  hx/2  (Ref  15 

108).  Since  the  difference  expression  error  is  0(hfc)  + 

2 

0(hx>,  obeying  the  stability  criteria  should  result  in  low 
error  with  neither  space  nor  time  error  dominating  the  total 


error  (Ref  15:108). 
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IV  Results  and  Compari sons 

The  results  of  MWR  and  finite  difference  solutions  to 
the  two  problems  are  compared  with  each  other  in  this  chap¬ 
ter.  The  results  are  also  compared  with  the  article  by 
Houstis  et.al.  which  compares  MWR  applied  on  finite  elements 
with  finite  differences  (Ref  6) . 

The  Houstis  et.al.  article  compares  collocation,  Gal- 
erkin,  least  squares,  and  finite  differences  for  linear 
second  order  elliptic  partial  differential  equations.  The 
MWR  methods  used  a  rectangular  grid  to  define  the  finite 
elements  and  Hermite  bicubic  polynomials  for  approximation 
of  the  solution.  Solution  of  equations  was  by  Gaussian 
elimination  by  profile  or  frontal  method.  The  article  con¬ 
clusions  state  collocation  was  more  efficient  than  finite 
differences  for  accuracy  of  one  to  four  significant  figures 
and  beyond.  The  measure  of  efficiency  was  accuracy  of 
solution  and  execution  time.  Finite  differences  and  collo¬ 
cation  started  with  equal  efficiency  at  low  accuracy  or, 
in  some  cases,  finite  differences  was  most  efficient  at  the 
lowest  accuracy.  As  accuracy  increased  collocation  became 
more  efficient  than  finite  differences  after  some  crossover 
point  at  one  to  four  significant  figures  of  accuracy. 
Accuracy  was  obtained  by  measuring  the  error  of  the  solution 
at  the  nodes  of  the  mesh  used  for  the  finite  element  and 
finite  difference  statement  of  the  problem  (Ref  10:323-334). 

The  comparison  of  Galerkin  and  least  squares  methods 
to  collocation  shows  that  collocation  is  always  faster  for 
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equal  accuracy  with  less  of  an  advantage  as  finer  grids  are 
used.  The  reason  for  the  speed  advantage  of  collocation  was 
attributed  to  its  narrower  band  matrix  than  Galerkin  or 
least  squares.  This  advantage  was  true  even  though  the  sym¬ 
metry  properties  of  the  Galerkin  and  least  squares  coeffi¬ 
cient  matrices  were  used  to  reduce  solution  time  by  half 
with  the  Cholesky  decomposition  for  band  matrices  (profile 
method) .  For  a  given  mesh  size  collocation  was  never  more 
accurate  than  the  other  methods.  The  main  conclusion  of 
the  article  is  that  collocation  is  the  best  method  of  the 
four  for  the  class  of  problems  examined  (Ref  6:335-337). 

Measurement  of  Accuracy 

The  measure  of  accuracy  used  for  comparisons  in  this 
study  is  the  maximum  absolute  error.  The  maximum  error 
encountered  in  all  points  tested  becomes  the  error  of  that 
solution.  Due  to  the  uniform  convergence  properties  of  MWR, 
the  maximum  error  is  expected  to  decline  with  more  terms  in 
the  trial  solution  as  finite  difference  error  declines  with 
more  mesh  points  (Ref  8:263-265,273).  The  fractional  error 
has  no  such  expectation.  Fractional  error  is  a  measure  of 
the  number  of  significant  figures  in  the  MWR  or  finite  dif¬ 
ference  solution.  Significant  figures  are  equal  to  -log 
(error/exact  solution)  where  error/exact  solution  is  frac¬ 
tional  error  (Ref  6:333).  Since  fractional  error  depends 
on  the  value  of  the  exact  solution  as  well  as  the  error  at  a 
point,  it  will  not  be  used  to  compare  efficiency  between 
methods. 
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Fractional  error  has  been  used  by  Houstis  et.al.  as  a 
characteristic  of  solution  comparisons.  After  finding  the 
crossover  point  for  finite  difference  and  collocation  a 
characteristic  number  of  digits  is  assigned  to  describe  the 
solutions  at  that  point.  At  the  crossover  point,  maximum 
error  and  execution  time  are  approximately  equal  with  collo¬ 
cation  obtaining  greater  accuracy  in  less  time  after  the 
crossover  point.  The  number  of  digits  assigned  to  the 
crossover  point  is  -log (maximum  error/maximum  solution  size). 
This  number  of  digits  is  clearly  an  optimistic  estimate  of 
the  number  of  significant  figures  in  the  solution  since  the 
fractional  error  of  the  maximum  error  is  computed  at  the 
largest  value  of  the  exact  solution  instead  of  where  the 
maximum  error  occurred  (Ref  6:333). 

The  absolute  error  must  be  measured  at  a  number  of 
points  in  the  domain  to  find  the  maximum  error.  Houstis  et. 
al.  measured  the  error  at  the  nodes  or  intersections  of  the 
finite  element  grid  (Ref  6:332).  These  are  the  mesh  points 
of  finite  difference  for  a  grid  of  the  same  size.  Since  no 
finite  element  mesh  is  used  for  the  problems  considered 
here,  the  finite  difference  points  will  be  used  for  compar¬ 
ison.  For  collocation  these  points  are  the  same  as  the 
collocation  points.  Then  for  an  N  term  expansion  of  the  MWR 
trial  solution,  the  error  will  be  measured  at  N  equally 
spaced  points,  initially.  A  more  extensive  investigation 
of  the  MWR  solutions  will  be  made  at  99  points  throughout 
the  domain  of  the  problem  to  check  the  error  obtained  at 

the  finite  difference  mesh  points. 
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Steady  State  Problem 


The  solution  to  the  steady  state  problem  is  shown  in 
Figure  1.  The  exact  solution  is  shown  in  comparison  to  the 
one  term  and  two  term  collocation  solutions.  The  conver¬ 
gence  of  the  MWR  solutions  to  the  true  solution  is  rapid 
with  successive  expansions  approximating  the  true  solution 
more  closely. 

In  order  to  compare  the  maximum  error  of  the  MWR  and 
finite  difference  solutions,  20  different  solutions  were 
computed  for  each  MWR  method  and  finite  differences.  Each 
MWR  method  used  solution  expansions  ranging  from  one  to  20 
terms.  The  finite  difference  solutions  ranged  from  one  to 
20  mesh  points.  The  results  of  maximum  error  and  execution 
time  measurements  are  tabulated  in  Table  I  and  Table  III 
through  Table  V  in  Appendix  A.  The  results  are  also  pre¬ 
sented  graphically  in  Figures  2  and  3. 

Figure  2  shows  the  maximum  absolute  error  obtained  by 
measurement  at  the  finite  difference  mesh  points  for  all 
four  methods  compared.  The  one  term  expansions  of  MWR  were 
compared  with  the  true  solution  at  one  point.  The  two  term 
expansions  were  compared  at  two  points  and  so  on.  The 
maximum  error  is  plotted  against  execution  time  for  100 
executions  of  each  solution.  Each  symbol  represents  a  dif¬ 
ferent  solution  with  a  different  number  of  terms  or  mesh 
points.  The  one  term  and  one  mesh  point  solutions  are  at 
the  top  of  each  line.  The  three  MWR  methods  performed 
similarly.  Galerkin  was  slightly  faster  than  least  squares 
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Maximum  Absolute  Error 


_  3 

and  collocation  m  the  maximum  error  range  around  10  to 

-  5  -  5 

10  with  least  squares  approaching  Galerkin  below  10 

Collocation  continues  slower  than  Galerkin  or  least  squares 

at  smaller  errors.  The  finite  difference  solutions  shown  in 

Figure  2  are  much  slower  to  a  given  maximum  error  than  the 

-4 

MWR  methods  for  errors  less  than  5x10  .  A  characteristic 

number  of  digits  accuracy  can  be  assigned  to  the  crossover 
point  for  this  problem  at  between  10  3  and  7x10  4  maximum 
error.  Since  the  maximum  solution  size  for  this  problem  is 
.071,  the  digits  of  accuracy  for  the  crossover  is  between 
1.8  and  2.0.  This  accuracy  is  a  very  optimistic  estimate 
of  the  accuracy  of  the  solution.  For  a  point  other  than  the 
crossover  different  accuracies  require  widely  varying  running 
times  by  different  methods.  At  maximum  error  of  7xl0-5  or 
3.0  digit  accuracy,  finite  differences  takes  twice  as  long 
as  the  MWR  methods. 

Figure  3  shows  maximum  error  plotted  against  execution 
time  for  100  executions  when  maximum  error  is  determined  by 
sampling  the  MWR  solutions  at  99  points.  The  absolute 
error  of  each  solution  was  measured  at  the  same  99  equally 
spaced  points  throughout  the  domain,  'the  finite  difference 
comparison  at  the  mesh  points  of  Figure  2  is  included  for 
reference.  Although  the  maximum  error  measured  for  MWR 
solutions  of  a  few  terms  is  a  few  to  several  times  higher 
for  99  points  than  for  the  finite  difference  points,  the 
difference  is  small  for  solutions  of  10  terms  or  more.  The 
efficiency  and  relative  efficiency  of  the  MWR  solutions  is 
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almost  the  same  for  measurement  at  the  99  points  as  the 
finite  difference  points. 

The  Houstis  et.al.  article  found  collocation  to  be 
faster  to  a  given  accuracy  than  Galerkin  or  least  squares 
for  MWR  applied  on  finite  elements  (Ref  6:335-337).  In  this 
formulation  of  MWR  on  the  whole  domain  and  for  this  problem, 
Galerkin  is  faster  than  either  least  squares  or  collocation. 
Also  least  squares  is  faster  than  collocation.  The  numeri¬ 
cal  integration  used  for  Galerkin  and  least  squares  in  MWR 
on  finite  elements  and  the  resulting  coefficient  matrix  may 
be  the  reasons  for  the  different  behavior.  As  indicated  in 
the  Houstis  article,  the  numerical  integration  for  MWR  on 
finite  elements  requires  evaluation  of  the  integrand  at  nine 
places  where  collocation  requires  evaluation  of  the  residual 
at  only  four  points  (Ref  6:336).  The  resulting  equations 
produce  coefficient  matrices  which  take  longer  to  solve  for 
Galerkin  and  least  squares  than  collocation  (Ref  6:336). 

MWR  on  the  whole  domain  as  used  in  this  study  has  neither 
of  these  two  hindrances.  The  integrals  are  done  analytically 
to  develop  recursion  relations  for  matrix  formation.  The 
resulting  matrices  have  all  non-zero  elements  and  are  all 
solved  by  the  same  method  so  solution  times  are  close  to 
the  same.  Galerkin  and  least  squares  were  observed  to  be 
more  accurate  for  a  given  mesh  size  than  collocation  by 
Houstis  et.al.  (Ref  6:336).  Apparently  the  greater  accuracy 
of  Galerkin  and  least  squares  results  in  less  error  for  a 
given  number  of  terms  in  the  expansion  than  collocation. 
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Then  the  very  close  solution  times  result  in  smaller  maximum 
error  for  Galerkin  than  for  collocation. 

Table  I  shows  some  of  the  values  presented  graphically 
in  Figure  2  and  Figure  3  for  collocation.  The  maximum  abso¬ 
lute  error  decreases  with  more  terms  in  the  trial  solution 
until  a  minimum  error  is  reached  at  11  or  12  terms.  After 
that  point  the  error  fluctuates  and  gradually  grows. 

Galerkin  and  least  squares  exhibit  similar  behavior  at 
eight  or  nine  terms  in  the  trial  solution  for  minimum  error. 

This  behavior  occurs  well  below  the  region  of  comparison  of 
Figures  2  and  3.  Galerkin  and  least  squares  error  are 

tabulated  in  Table  III  and  Table  IV  of  Appendix  A.  All 

-12 

three  solutions  attained  minimum  errors  of  10  or  less. 

There  are  three  possible  explanations  for  the  behavior 
of  maximum  error.  The  first  explanation  of  convergence  to 
some  solution  other  than  the  true  solution  is  not  valid. 

If  the  three  methods  converged  to  the  wrong  solution,  maxi¬ 
mum  error  would  decrease  to  some  minimum  value  and  not  grow. 

A  second  explanation  is  low  accuracy  of  the  Gaussian  elimi¬ 
nation  solution  of  the  simultaneous  equations.  A  successive 
over-relaxation  iterative  method  was  used  to  improve  the 
Gaussian  elimination  solution  (Ref  15:126-129).  A  relative  con¬ 
vergence  test  of  10  ^  was  used  to  test  for  convergence  of 
the  iterations.  For  collocation  the  successive  over-relax¬ 
ation  gave  greater  error  for  some  trial  solutions  and  failed 
to  yield  a  solution  for  others.  For  Galerkin  and  least 
squares,  however,  the  method  worked.  Several  solution  con¬ 
stants  changed  value  in  the  10th  to  11th  significant  figures. 
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TABLE  I 


Error  and  Time  for  Collocation  on  the 
Steady  State  Problem 


Maximum  Absolute  Error 


Number  of  Terms  At 

in  the  Collociition  At  99 

Trial  Solution  Points  Points 


Time  for 
100  Executions 
_ (seconds ) 


1 

2 

3. 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 


1. 68E-3 

7.09E-4 

4. 87E-5 

3. 23E-6 

1.96E-7 

7.73E-9 

4.08E-10 

1.16E-11 

5.35E-13 

1.46E-14 

3.77E-15 

5.77E-15 

8.44E-15 

4.22E-15 

6 . 22E-15 

6. 22E-15 

4.85E-14 

2.29E-14 

5.77E-14 

1. 38E-13 


9.69E-3 
8 . 06E-4 
7.06E-5 
3. 38E-6 
2. 38E-7 
7. 93E-9 
4.65E-10 
1.18E-11 
5.90E-13 
1. 88E-14 
6 . 77E-15 
6.55E-15 
9.77E-15 
5. 33E-15 
7.11E-15 
6. 66E-15 
4.88E-14 
2. 49E-14 
5. 99E-14 
1. 38E-13 


.18 
.  21 
.25 
.30 
.  37 
.45 
.  53 
.62 
.71 
.85 
.98 
1.11 
1.25 
1.42 
1.62 
1.81 
2.01 
2.23 
2.45 
2.65 


Note:  E-x  means  10-X. 
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The  maximum  errors  of  the  solutions  in  Tables  III  and  IV 
remained  the  same  to  three  significant  figures  yielding 
Tables  III  and  IV  again.  Since  the  growth  in  the  error  was 
up  to  two  orders  of  magnitude,  poor  accuracy  of  the  Gaussian 
elimination  does  not  explain  the  behavior  of  maximum  error. 
The  last  explanation  of  the  maximum  error  behavior  is  the 
trial  functions  themselves.  Babuska  et.al.  find  these 
trial  functions  numerical  unstable  for  the  MWR  in  their 
study  of  stability  in  optimal  trial  functions  (Ref  14:241- 
245) .  By  studying  the  solution  of  the  simultaneous  equa¬ 
tions  based  on  these  trial  functions,  they  conclude  they  are 
unstable  for  seven  or  more  terms  in  the  trial  solution.  The 
solutions  with  more  than  seven  terms  have  errors  which 
increase,  not  decrease.  Babuska  et.al.  conclude  these  trial 
functions  should  not  be  used  for  computer  solutions. 

Transient  Problem 

The  transient  problem  was  solved  by  the  reduction  to 
ordinary  differential  equations  method  of  MWR  and  by  an 
explicit  finite  difference  method.  The  MW7R  on  the  whole 
domain  produces  a  series  solution  that  may  be  evaluated  for 
any  time  or  position  in  the  domain.  The  explicit  finite 
difference  method,  however,  produces  the  solution  for  a 
succession  of  times  up  to  the  final  time  considered.  For 
the  comparison  of  the  methods  in  this  study,  only  the  maxi¬ 
mum  absolute  error  at  specified  dimensionless  times  will  be 
used . 
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Several  points  within  the  space  domain  were  used  to 
compare  an  approximate  solution  and  the  exact  solution.  The 
finite  difference  solutions  were  found  for  10  to  100  mesh 
points  in  increments  of  10.  The  finite  difference  mesh 
points  were  used  to  determine  the  maximum  absolute  error  of 
that  solution.  The  MWR  solutions  used  successively  more 
terms  and  solutions  with  1  to  20  terms  were  calculated.  The 
MWR  solutions  were  compared  at  the  collocation  points  for 
the  steady  state  problem.  Since  most  of  the  finite  differ¬ 
ence  solutions  in  the  transient  problem  were  compared  at 
more  than  the  maximum  number  of  collocation  points,  all  the 
MWR  solutions  will  be  compared  at  99  equally  spaced  points 
for  a  comparable  search  of  the  domain. 

The  solutions  were  compared  at  three  times.  Figure  4 
shows  the  exact  solutions  at  t  =  .05,  .10,  and  .15.  Note 
that  the  solution  for  t=.05  is  close  to  zero  for  a  larger 
fraction  of  the  domain  than  the  others  and  represents  the 
earliest  time  response  of  the  solution.  The  solution  for 
t=.15  represents  the  latest  time  response  when  the  solution 
is  close  to  the  steady  state  solution. 

The  maximum  absolute  error  for  the  three  solutions  is 
plotted  as  a  function  of  solution  execution  time  in  Figures 
5  through  7.  All  solutions  were  executed  100  times  for  a 
more  accurate  measurement  of  the  time.  Only  collocation 
and  Galerkin  MWR  methods  were  us  '  lalerkin  was  always 
faster  to  a  given  accuracy  than  collocation  as  in  the 
steady  state  problem.  Both  the  Galerkin  and  collocation 
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Fig  4.  Transient  Prob; 


□  -  Finite  Differences 
Q-  Galerkin 
A-  Collocation 


Fig  5.  Transient  Error  for  Time  =  .05 
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Maximum  Absolute  Error 
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methods  were  faster  to  a  given  accuracy  than  explicit  finite 
differences  after  a  crossover  point  where  MWR  and  finite 
differences  achieved  equal  accuracy  in  equal  time.  The 
crossovers  occurred  at  accuracies  of  2.5  to  3.2  digits  based 
on  the  maximum  absolute  error  and  maximum  solution  size  of 
1.0.  The  crossovers  represent  very  optimistic  estimates  of 
accuracy.  All  methods  showed  improved  accuracy  for  t=.10 
and  t=.15  over  t=.05.  Collocation  improved  the  most  and 
finite  differences  the  least.  The  crossover  point  for 
collocation  and  finite  differences  changed  from  3.2  digits 
at  t= . 05  to  2.6  digits  at  t=.15.  The  Galerkin  and  finite 
differences  crossover  remained  more  nearly  constant  by 
changing  from  2.7  digits  at  t=.05  to  2.5  digits  at  t=.15. 
Collocation  lost  more  of  its  advantage  over  finite  differ¬ 
ences  at  the  short  solution  time  t=.05  than  did  Galerkin. 

An  explanation  for  the  better  showing  of  Galerkin  over 
collocation  is  the  complex  eigenvalues  found  for  colloca¬ 
tion.  The  Galerkin  matrices  in  the  eigenvalue  determina¬ 
tion  part  of  the  solution  process  were  symmetric,  while  the 
collocation  matrices  were  not.  For  collocation  complex 
eigenvalues  were  first  found  for  the  five  term  expansion  of 
the  trial  solution  and  expansions  for  more  terms  also  had 
complex  eigenvalues.  Since  only ■ the  real  part  of  a  complex 
eigenvalue  was  used  as  an  approximation  in  the  solution, 
some  accuracy  was  lost.  Accuracy  was  also  lost  for  the 
Galerkin  method  due  to  negative  real  eigenvalues.  These 
negative  eigenvalues  produced  positive  exponential  functions 
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of  time  which  grew  instead  of  decayed  with  time.  Since  the 
negative  eigenvalues  produced  functions  which  were  not  a 
physically  real,  decaying  transient  solution,  those  eigen¬ 
values  were  set  to  zero.  The  negative  eigenvalues  first 
occurred  for  the  11  term  Galerkin  solution  and  for  all 
solutions  with  more  terms.  Thus,  the  effect  on  accuracy 
of  the  negative  eigenvalues  for  Galerkin  occurred  at  finer 
approximations  than  the  complex  eigenvalues  for  collocation. 

An  explanation  for  the  poorer  showing  of  both  MWR 
methods  at  the  t=.05  time  is  the  shape  of  the  true  solution 
at  that  time.  The  true  solution  there  is  close  to  zero  for 
a  greater  fraction  of  the  domain  than  the  other  two  times. 
Since  the  true  solution  is  close  to  zero  for  that  longer 
fraction  of  the  domain,  it  is  harder  to  approximate  by  a 
polynomial.  More  terms  and  more  solution  time  are  required 
for  a  given  accuracy  or  maximum  absolute  error. 

Table  II  shows  the  maximum  absolute  error  and  execu¬ 
tion  times  for  collocation  solutions  at  the  three  dimension¬ 
less  times.  Note  that  the  maximum  absolute  error  stops 
decreasing  and  starts  growing  for  solutions  of  13  or  14 
terms.  The  error  for  Galerkin  also  grows  starting  at  13 
terms  in  the  solution  as  shown  in  Table  VI  of  Appendix  A. 
Since  the  same  space  trial  functions  were  used  for  the  solu¬ 
tion  expansions  in  the  steady  state  and  transient  problems 
and  since  similar  growth  in  error  occurred  in  the  steady 
state  problem,  this  growth  in  error  can  be  attributed  at 
least  partially  to  the  numerically  unstable  trial  functions 


54 


TABLE  II 


Error  and  Time  for  Collocation  in  the 
Transient  Problem 


Maximum  Absolute  Error 


Number  of  Terms 
in  the 

Trial  Solution  T=.05  T=.10  T=.15 


Time  for 
100  Executions 
(seconds) 


1 

7.69E-2 

1. 39E-2 

1.21E-2 

.90 

2 

2.76E-2 

1.75E-2 

1. 15E-2 

1.34 

3 

2.17E-2 

6.15E-3 

1. 92E-3 

2.11 

4 

6.66E-3 

1. 67E-3 

7. 68E-4 

3.07 

5 

2.79E-3 

5.55E-4 

1.63E-4 

5.22 

6 

2.24E-3 

7.37E-5 

3. 89E-5 

8.13 

7 

4. 57E-4 

3.69E-5 

8.09E-6 

10.69 

8 

1.95E-4 

5. 33E-6 

8.71E-7 

14.76 

9 

3.88E-5 

1.71E-6 

3.01E-7 

19.46 

10 

1.51E-5 

3.87E-7 

1. 88E-8 

26.89 

11 

4.91E-6 

8.98E-8 

8.86E-9 

34.20 

12 

1.24E-6 

2.07E-8 

7.07E-10 

38.95 

13 

4. 52E-7 

4.13E-9 

8.43E-10 

47.26 

14 

7.82E-8 

1.40E-8 

8.83E-9 

60.19 

15 

1.16E-7 

6.38E-8 

3.87E-8 

72.17 

16 

5.35E-7 

2.67E-7 

1.60E-7 

82.23 

17 

2.48E-6 

1.53E-6 

9.33E-7 

— 

18 

9.07E-6 

2 . 06E- 6 

8. 51E-7 

— 

19 

3.74E-5 

2.36E-5 

1.44E-5 

— 

20 

2.35E-2 

4. 34E-3 

8. 00E-4 

—  —  — 

Note : 

E-x  means 

10~X. 

(Ref  14:241-245).  No  iterative  improvement  of  the  final 
solution  was  attempted,  however. 

An  implicit  finite  difference  method  was  used  to  see 

if  the  solution  time  could  be  shortened  from  the  explicit 

2 

method.  Using  a  time  step,  hfc  =  h^,  where  hfc  is  the  time 
step  size  and  where  h^  is  the  space  mesh  spacing  to  minimize 
the  finite  difference  error,  0(ht)  +  O(h^)  (Ref  14:108). 

The  resulting  execution  times  were  several  to  100  times 
greater  than  for  the  explicit  solutions  with  slightly  higher 
maximum  absolute  error.  The  largest  increase  in  solution 
times  occurred  for  the  finest  mesh  spacing  where  the  largest 
full  matrices  were  solved  for  a  new  temperature  more  times 
than  the  smaller  matrices  for  larger  mesh  sizes. 

No  direct  comparison  with  the  Houstis  et.al.  article 
is  possible  for  the  transient  problem  since  that  article 
treated  only  steady  state  problems  on  finite  elements 
(Ref  6) .  It  is  interesting  to  note  that  the  Galerkin 
method  proved  faster  to  a  given  accuracy  than  collocation 
for  both  problems  in  this  study  while  collocation  was  fastest 
for  the  MWR  on  finite  elements  of  the  Houstis  et.al.  article. 
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V  Conclusions 


In  this  chapter  conclusions  of  the  investigation  of 
MWR  on  the  whole  domain  and  finite  differences  for  solution 
of  the  heat  equation  will  be  presented.  Also,  recommenda¬ 
tions  for  the  improvement  of  the  MWR  methods  used  in  this 
study  and  application  of  finite  element  methods  to  heat 
problems  of  axial  symmetry  will  be  presented. 

Two,  one  dimensional  problems  of  heat  transfer  were 
solved  by  the  MWR  on  the  whole  domain  of  interest  and  by 
finite  differences.  One  problem  represented  steady  state 
heat  transfer  and  the  other  transient  heat  transfer.  Both 
problems  used  Dirichlet  boundary  conditions  where  the  value 
of  temperature  was  specified  on  the  boundary.  MWR  methods 
for  the  steady  state  problem  were  collocation,  Galerkin, 
and  least  squares.  The  least  squares  method  was  not  used 
for  the  transient  problem  since  it  complicated  the  solution 
process.  The  finite  difference  formulation  was  implicit 
for  the  steady  state  problem  and  explicit  for  the  transient 
problem.  The  integrals  of  the  Galerkin  and  least  squares 
methods  were  evaluated  analytically  to  obtain  recursion 
relations  for  MWR  solution.  All  solutions  were  compared  with 
the  true  solution  to  obtain  maximum  absolute  error. 

Conclusions  and  Limitations 

The  Galerkin  method  was  fastest  to  a  given  accuracy 
or  maximum  absolute  error  of  the  methods  evaluated.  The 
crossover  point  for  Galerkin  and  finite  differences  where 
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both  had  equal  accuracy  in  equal  time  occurred  at  1.8  to  2.0 
digits  of  accuracy  for  the  steady  state  problem  and  2.5  to 
2.7  digits  for  the  times  evaluated  in  the  transient  problem. 
The  Galerkin  method  had  greater  accuracy  than  finite  dif¬ 
ferences  in  equal  time  after  the  crossover  point.  The 
digits  of  accuracy  are  an  optimistic  estimate  of  the  number 
of  significant  figures  in  the  approximate  solution  based  on 
the  maximum  absolute  error  occurring  in  the  whole  domain 
and  the  maximum  value  of  the  true  solution.  The  digits  of 
accuracy  indicate  that  finite  differences  is  only  faster  for 
rather  crude  accuracy  and  Galerkin  and  the  other  MWR 
methods  are  faster  for  greater  accuracy.  The  speed  margin 
is  sizeable.  The  time  required  to  increase  one  digit  of 
accuracy  beyond  the  crossover  point  for  the  steady  state 
problem  is  half  as  long  for  Galerkin  as  for  finite  differ¬ 
ences.  The  Galerkin  method  is  the  best  of  methods  evalu¬ 
ated  to  use  for  problems  of  this  type  when  accuracy  greater 
than  two  or  three  digits  is  required. 

Some  limitations  to  the  conclusion  of  the  previous 
paragraph  must  be  noted.  First  only  Dirichlet  boundary 
conditions  were  used  for  the  two  problems  evaluated.  Other 
boundary  conditions  may  be  evaluated  by  finding  trial  func¬ 
tions  which  satisfy  them  or  by  minimizing  the  boundary 
residual  along  with  the  differential  equation  as  shown  in 
Eq.  (10)  of  Chapter  II.  However,  the  conclusion  of  the 
last  paragraph  cannot  be  extended  to  other  boundary  condi¬ 
tions  without  some  evaluation  of  problems  with  the  other 
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bo  u ' :  ary  conditions.  Another  limitation  is  constant  thermal 
conductivity.  Constant  thermal  conductivity  throughout  the 
whole  domain  was  assumed  to  reduce  the  heat  equation  to  a 
dimensionless  form.  Problems  involving  varying  thermal 
conductivity  should  be  evaluated  to  see  if  Galerkin  is  the 
best  method.  A  third  limitation  is  the  one  dimensional 
cartesian  coordinate  geometry  of  both  problems  evaluated. 
Before  assuming  that  Galerkin  is  always  fastest  for  axial 
symmetry  or  two  dimensional  problems  on  the  whole  domain 
problems  in  those  geometries  should  be  evaluated.  The  last 
limitation  is  the  method  of  integral  evaluation  used  for 
the  Galerkin  and  least  squares  methods.  If  a  numerical 
rather  than  an  analytical  integration  were  used,  the  solu¬ 
tion  times  and,  therefore,  the  speed  advantage  of  those  two 
methods  could  change. 

A  major  drawback  of  the  Galerkin  and  other  MWR  methods 
on  the  whole  domain  in  the  problems  evaluated  is  the  poly¬ 
nomial  trial  functions.  These  trial  functions  were  seen  to 
be  numerically  unstable  for  solutions  with  10  or  more  terms. 
The  maximum  absolute  error  increased  not  decreased  for 
solutions  with  more  than  10  terms  or  so.  If  heat  generation 
terms  complicated  the  shape  of  the  temperature  distribution, 
more  terms  could  be  needed  in  the  MWR  solution  to  approxi¬ 
mate  it.  If  more  than  10  terms  were  needed,  MWR  might  be 
limited  to  low  accuracy  and  be  slower  than  finite  differ¬ 
ences.  A  new  approach  is  needed  to  avoid  the  drawback  of 
the  linearly  independent  but  numerically  unstable  polynomial 
trial  functions. 
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Recommendations 


The  orthogonal  collocation  method  could  solve  the  prob¬ 
lem  of  numerical  instability  (Ref  3:357).  it  is  recommended 
for  steady  state  problems.  Orthogonal  collocation  uses 
trial  functions  composed  of  orthogonal  polynomials  and 
collocation  points  defined  by  the  roots  of  the  highest 
order  polynomial.  The  method  can  be  made  to  fit  other  bound¬ 
ary  condition  than  Dirichlet  and  is  valid  for  planar,  cylin¬ 
drical  or  spherical  geometry.  The  method  can  solve  problems 
in  terms  of  the  values  of  the  function  approximated  at  the 
collocation  points.  This  ability  allows  great  flexibility 
as  well  as  computational  ease  in  problem  solution.  Ortho¬ 
gonal  collocation  has  been  shown  to  be  as  fast  as  the  other 
MWR  methods  for  a  given  accuracy  (Ref  3:97,100). 

As  an  example  consider  a  one-dimensional  cylindrical 
symmetry  steady  state  heat  transfer  problem  where  tempera¬ 
ture  varies  only  with  the  radius: 

V2U(r)  +  g (r) U (r)  =  f (r)  (85) 

0  <  r  <  1 

where  U(r)  represents  dimensionless  temperature  and  where 
g(r)  and  f (r)  are  functions  defining  heat  generation.  The 
problem  can  he  formulated  as  symmetric  in  r: 

V2U (r )  +  g (r2) U (r)  =  f(r2)  (86) 

Then  the  trial  solution  for  such  a  problem  can  be  formulated 
for  Dirichlet  boundary  conditions  at  r=l  as 
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U (r )  =  U(l)  +  (1-r2)  X  aiP._1(r2)  (87) 

i=l 

2 

where  is  an  undetermined  constant  and  P^_^ (r  )  is  an 

2 

orthogonal  polynomial.  The  P^_^  (r  )  are  defined  starting 

2 

with  a  power  series  in  r  as 

1 

w  (r2)  P  j  (r2)  Pi  (r2)  r  dr  =  (88) 

2 

for  j=l , 2 , . . . , i-1  where  w(r  )  is  a  weight  function,  6^^  is 

the  Kronecker  delta  function,  and  some  constant.  Thus  each 

polynomial  in  the  trial  function  of  Eq.  (87)  is  made  ortho- 

2 

gonal  to  the  others  with  the  weight  function,  w(r  ).  The 

roots  of  the  polynomials  may  be  found  for  a  given  weight 

2 

function.  Weight  function  (1-r  )  or  1  may  be  used  for  this 
problem.  The  polynomial  roots  for  these  weight  functions 
have  been  tabulated  (Ref  3:99,101-103). 

For  solution  the  trial  function  must  be  converted  to 
an  ordinary  polynomial  and  matrices  found  to  relate  the 

values  of  temperature  to  the  derivatives  in  the  problem. 

2  2 
Since  PN_^(r  )  is  a  polynomial  of  degree  N-l  in  r  ,  then 

the  trial  function  can  be  represented  as  a  polynomial  of 
2 

degree  N  in  r  : 

N+1  2i-2 

U(r)  =  X  d.r  (89) 

i=l  1 

where  d^  are  new  constants.  Evaluating  the  derivative  and 
Laplacian  of  temperature  at  the  collocation  points  plus 
r=l  for  the  N+1  collocation  point  yields  matrices  to  relate 
temperature  and  derivatives  to  the  constants  d^: 
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N+l  _ 

U  (r  .  )  =  >:  x?W  d. 

3  i=l  J  1 


(90) 


dU 

dr 


r  . 
3 


mi  dr2i-2 
i=l 


dr 


d. 

i 


(91) 


r  . 
3 


7 

V  U 


r  . 
3 


N+l  _  - .  - 
I  V2r21"2 
i=l 


(92) 


r  . 
3 


In  matrix  form  Eq.  (90)  is  U  =  Q  d,  Eq.  (91)  is  ~  =  C  d, 

2 

and  Eq.  (92)  is  7  U  =  D  d.  The  derivative  and  Laplacian 
can  now  be  formulated  in  terms  of  the  temperature  at  the 
collocation  points: 


dU 

dr 


=  C  Q  1  U 


(93) 


V2U  =  D  Q-1  U 


(94) 


The  expressions  of  Eq.  (93)  and  (94)  may  be  substituted 
directly  into  a  differential  equation  such  as  Eq.  (86)  to 
produce  a  matrix  problem  for  temperature  at  the  collocation 
points.  Then  the  temperature  at  the  collocation  points  may 
be  used  to  solve  for  the  constants,  d.,  in  the  trial  solu- 

l 

tion.  An  additional  condition  needed  for  solution  of  Eq. 

(86)  is  the  temperature  at  r=l  or  U  ( r^ )  =  U(l).  Other 
boundary  conditions  can  be  solved  for  u(rN+^)  an(3  included 
in  the  solution  for  temperature  (Ref  3:100-101). 

Another  method  which  should  be  evaluated  in  future 
research  is  the  MWR  on  finite  elements  for  problems  of  axial 
symmetry.  MWR  on  finite  elements  has  been  shown  to  be 
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faster  to  a  given  accuracy  than  finite  differences  for 
steady  state  problems  similar  to  the  heat  equation  in  two 
dimensions  (Ref  5:336).  Finite  elements  methods  use  piece- 
wise  continuous  trial  functions  of  a  fixed  number  of  terms. 
Accuracy  is  increased  by  increasing  the  number  of  elements 
covering  the  domain  of  the  problem  and  not  the  number  of 
terms  in  the  trial  function. 

Problems  of  axial  symmetry  use  volume  finite  elements 
defined  in  the  two  variables  r  and  z.  Axial  symmetry  is 
described  by  a  cylindrical  coordinate  system  where  tempera¬ 
ture  or  boundaries  do  not  vary  with  the  angular  coordinate 
of  the  cylindrical  coordinate  system.  The  domain  of  the 
problem  is  then  divided  into  finite  element  rings  with 
cross  section  and  trial  functions  described  by  r  and  z 
coordinates  of  the  cylindrical  coordinate  system.  For  inte¬ 
grations  involved  in  the  MWR  the  volume  of  the  finite  ele¬ 
ment  ring  must  be  used.  Various  shape  elements  may  be  used 
with  the  appropriate  trial  functions  to  element  shape  and 
order  of  the  polynomial  trial  function  used.  The  trial 
solution  is  expressed  as  in  orthogonal  collocation  as  a 
function  of  parameters  at  specific  nodes  or  points  on  the 
element : 

U  =  I  N.  a.  (95) 

where  U  is  the  function  approximated  such  as  temperature, 
are  the  trial  functions  called  shape  functions  in  finite 
elements  dependent  on  the  shape  and  order  of  polynomial 
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approximation,  and  where  are  the  nodes  determined  by  ele¬ 
ment  shape  and  shape  functions  (Ref  Id : 119-120, 148-177) . 

The  MWR  statement  of  the  steady  state  heat  transfer 
problem  minimizes  the  combined  residual  of  the  differential 
equation  and  boundary  conditions  of  second  and  third  kind 
to  produce  a  matrix  statement  for  solution  of  the  nodal 
parameters,  a^ .  When  the  matrix  equations  are  solved  the 
solution  is  complete-  The  matrix  equations  are  of  the  form 
H  a  +  £  =  0  where  H  is  obtained  by  numerical  integration 
of  the  integrals  involving  a  and  where  f  is  obtained  from 
integrals  not  involving  a  (Ref  13:424-426).  The  method  of 
combining  contributions  to  the  H  matrix  from  each  finite 
element  of  the  domain  results  in  a  banded,  sparse  matrix 
which  is  faster  to  solve  than  the  full  matrices  of  MWR  on 
the  whole  domain  (Ref  13:14-15). 

The  transient  heat  transfer  equation  may  be  treated  in 
the  same  manner  by  MWR  to  produce  ordinary  differential 
equations  in  time.  The  matrix  formation  is  then  of  the 
form  Ca+Ha+f=  0  which  may  be  solved  by  methods  simi¬ 
lar  to  those  used  for  the  transient  problem  in  this  study. 

A  step  by  step  recurrence  calculation  for  a^  similar  to 
finite  difference  methods  for  the  transient  problem  is  also 
possible  and  more  general  (Ref  13:569).  In  the  step  by  step 
method  shape  functions  are  chosen  to  describe  the  variation 

of  a  from  the  beginning  to  the  end  of  a  time  element, 

n+1 

a  =  E  N.  a..  The  MWR  is  applied  to  the  whole  matrix 

i-n  1  1 

formulation  of  the  problem  with  the  time  derivatives  acting 
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on  the  shape  functions.  The  resulting  equation  is  then 
reduced  to  yield  a  recurrence  relationship  for  a_n+^  as  a 
function  of  a^.  Thus  the  solution  can  be  produced  in  a 
step  by  step  basis  as  for  finite  differences  (Ref  13:570-572). 
Evaluation  of  the  finite  element  methods  described  here 
should  provide  a  better  and  more  general  test  of  MWR  for 
heat  transfer  problems  of  axial  symmetry. 
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Appendix  A:  Error  and  Solution  Time  Tables 

This  appendix  contains  tables  summarizing  the  maximum 
absolute  error  and  solution  time  for  MWR  and  finite  differ¬ 
ence  solutions  of  the  steady  state  and  transient  heat 
transport  problems. 


68 


TABLE  III 


Error  and  Time  for  Galerkin  on  the 
Steady  State  Problem 


Maximum  Absolute  Error 


Number  of  Terms 
in  the 

Trial  Solution 


At  Finite 
Difference 
Mesh  Points 


Time  for 

At  99  100  Executions 

Points  (seconds) 


1 

3. 03E-4 

8.35E-3 

.18 

2 

1.49E-4 

3.04E-4 

.20 

3 

1. 87E-5 

2. 57E-5 

.23 

4 

5. 22E-7 

5. 88E-7 

.27 

5 

3.30E-8 

3.74E-8 

.  30 

6 

6 . 19E-10 

6. 40E-10 

.36 

7 

2. 89E-11 

3.21E-11 

.42 

8 

8.10E-13 

8.14E-13 

.46 

9 

5.43E-12 

5. 52E-12 

.52 

10 

2.29E-11 

2.50E-11 

.58 

11 

3.80E-11 

3.80E-11 

.67 

12 

4.60E-11 

4.75E-11 

.75 

13 

3.72E-11 

3. 78E-11 

.81 

14 

1.71E-10 

1.74E-10 

.91 

15 

1.69E-10 

1.81E-10 

1.00 

16 

6.45E-11 

6.68E-11 

1.12 

17 

6.13E-11 

6.34E-11 

1.25 

18 

1.36E-10 

1.41E-10 

1.37 

19 

9.86E-11 

1.09E-10 

1.47 

20 

1.03E-10 

1.17E-10 

1.62 

Note:  E-x  means  10 


-x 
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TABLE  IV 


Error  and  Time  for  Least  Squares  on  the 
Steady  State  Problem 


Maximum  Absolute  Error 


Number  of  Terms  At  Finite  Time  for 

in  the  Difference  At  99  100  Executions 

Trial  Solution  Mesh  Points  Points  (seconds) 


1 

1. 68E-3 

9 . 06E-3 

.19 

2 

1.42E-3 

1.68E-3 

.20 

3 

7.80E-5 

7. 90E-5 

.24 

4 

1.33E-6 

1.85E-6 

.  28 

5 

9.60E-8 

9. 87E-8 

.35 

6 

1.33E-9 

1.68E-9 

.39 

7 

7. 53E-11 

7.84E-11 

.49 

8 

9.94E-13 

1.06E-12 

.55 

9 

4.56E-13 

4.84E-13 

.65 

10 

8.58E-13 

8.77E-13 

.75 

11 

3.96E-12 

3.98E-12 

.86 

12 

2.71E-11 

2.73E-11 

.97 

13 

1.01E-11 

1.03E-11 

1.09 

14 

2.76E-11 

3.14E-11 

1.23 

15 

1.29E-11 

1.31E-11 

1.35 

16 

3.94E-11 

4.12E-11 

1.53 

17 

7.41E-11 

7.48E-11 

1.71 

18 

2.09E-10 

2.40E-10 

1.84 

19 

6.44E-11 

6.73E-11 

2.05 

20 

3.65E-11 

3.77E-11 

2.22 

Note:  E-x  means  10 
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TABLE  V 


Error  and  Time  for  Finite  Difference  Method 
on  the  Steady  State  Problem 


Maximum  Absolute  Error 

Number  of  Terms  At  Finite  Time  for 

in  the  Difference  100  Executions 

Trial  Solution  Mesh  Points  (seconds ) 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 


1. 68E-3 
7.09E-4 
4 . 09E-4 
2. 64E-4 
1. 81E-4 
1.35E-4 
1.02E-4 
8.17E-5 
6. 58E-5 
5.46E-5 
4. 58E-5 
3.90E-5 
3.37E-5 
2.92E-5 
2.58E-5 
2. 28E-5 
2.04E-5 
1.83E-5 
1.65E-5 
1.50E-5 


.26 
.23 
.26 
.  32 
.36 
.43 
.47 
.54 
.60 
.  66 
.76 
.84 
.89 
.98 
1.09 
1.20 
1.32 
1.41 
1.53 
1.67 


Note:  E-x  means  10 


TABLE  VI 


Error  and  Time  for  Galerkin  in  the 
Transient  Problem 


Maximum  Absolute  Error 


Number  of  Terms 
in  the 

Trial  Solution  T= . 0 5  T= . 10  T~ . 15 


Time  for 
100  Executions 
(seconds) 


1 

5. 39E-2 

1.45E-2 

6. 18E-3 

.91 

2 

1.92E-2 

9.86E-3 

5.65E-3 

1.33 

3 

6. 60E-3 

1. 56E-3 

3. 76E-4 

1.99 

4 

1.19E-3 

2. 25E-4 

1. 34E-4 

2.96 

5 

4.52E-4 

6.98E-5 

1.08E-5 

4.57 

6 

1.09E-5 

3. 00E-6 

1. 53E-6 

7.03 

7 

1.43E-6 

2.41E-6 

3. 38E-7 

9.51 

8 

7.01E-7 

9. 92E-8 

1. 32E-8 

12.18 

9 

4.82E-7 

5. 53E-8 

7. 65E-9 

16.01 

10 

3.84E-7 

6. 81E-9 

2. 24E-9 

20.32 

11 

5. 76E-8 

5.45E-9 

1. 98E-9 

24.80 

12 

1.70E-7 

2. 07E-9 

4. 44E-10 

31.10 

13 

3.95E-8 

1.35E-9 

2.21E-10 

38.47 

14 

8.13E-7 

4.04E-8 

2 . 80E-9 

42.69 

15 

1.90E-7 

1. 26E-8 

9. 78E-10 

50.99 

16 

4.18E-7 

2. 66E-7 

1. 68E-9 

59.39 

17 

6. 27E-7 

1.71E-7 

4.20E-8 

68.25 

18 

5.75E-7 

1.41E-7 

4.62E-8 

77.08 

19 

5. 94E-7 

1.41E-7 

3.00E-8 

— 

20 

1.11E-5 

5. 83E-6 

3. 09E-6 

Note : 

E-x  means 

10_X. 
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TABLE  VII 


Error  and  Time  for  the  Explicit  Finite  Difference 
Method  in  the  Transient  Problem 

T=.  05 

Time  for 

Number  of  Maximum  100  Executions 

Mesh  Points  Absolute  Error  (seconds) 

10 

1.10E-2 

.69 

20 

2.74E-3 

1.77 

30 

1.28E-3 

3.96 

40 

7. 30E-4 

7.66 

50 

4.72E-4 

13.45 

60 

3. 31E-4 

21.59 

70 

2. 44E-4 

32.81 

80 

1. 88E-4 

47.00 

90 

1.49E-4 

65.45 

100 

1.21E-4 

88.21 

Note : 

E-x  means  10  x.  T 

is  the 

dimensionless  time 

variable. 

Error  and  Time  for 
Method  in 

Number  of 

Mesh  Points  Ab 

TABLE  VIII 

the  Explicit 
the  Transient 
T=.  10 

Maximum 
solute  Error 

Finite  Difference 
Problem 

Time  for 

100  Executions 
(seconds) 

10 

5.13E-3 

.85 

20 

1.41E-3 

2.69 

30 

6. 53E-4 

6.73 

40 

3. 73E-4 

13.71 

50 

2. 4  IE- 4 

24.62 

60 

1. 69E-4 

40.46 

70 

1. 25E-4 

62.44 

80 

9. 58E-5 

91.17 

90 

7. 59E-5 

127.52 

100 

6.16E-5 

— 

Note:  E-x  means  10  X.  T  is  the 

dimensionless  time  variable. 
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TABLE  IX 


Error  and  Time  for  the  Explicit  Finite  Difference 
Method  in  the  Transient  Problem 
T=.  15 


Number  of 
Mesh  Points 

10 

20 

30 

40 

50 

60 

70 

80 

90 

100 

Note : 


Maximum 
Absolute  Error 

3. 65E-3 
1. 02E-3 
4. 69E-4 
2.68E-4 
1.73E-4 
1.21E-4 
8. 96E-5 
6. 88E-5 
5. 45E-5 
4 . 43E-5 


Time  for 
100  Executions 
( seconds ) 


1.02 

3.55 

9.25 

19.57 

35.56 

58.99 

91.43 

135.13 


E-x  means  10  x.  T  is  the 
dimensionless  time  variable. 


Appendix  D:  Computer  Programs 


This  appendix  contains  the  computer  programs  for  the 
steady  state  and  transient  heat  transfer  problems.  Doth 
the  MWR  and  finite  difference  solutions  are  included. 
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Steady  State  Computer  Program 


PROGRAM  EASY!  ( INPUT, GUTPUT) 

DIP! '.SIGN  CCE1'  (  20 , 20 ) , RS ( 20)  >  «K AREA ( 20 )  ,  F D  ( 20 ) 
I D I  Y  =  2  0 

:o  i  icri t= 1 , 3 

LD  :  N  3  1 ,20 
N 1  -  \  + 1 

CAL -  KKRSK ( CGEr  -  RS ,  wK AREA ,N,  IDIM,  ICRIT) 

CALL  Z  L  M  P  M  a  (  R  S  ,  N  ,  I  D  I  M ) 

CAL-  FDSM(CCEF,FD,XKAREA,N1 , IDIM) 

CA_L  CGMPFD(FD.Nl.IDIM) 

1  CONTINUE 
END 

SUBROUTINE  MaRSM ( CCEE ,  R5,kKAR’EA,N, IDIM, ICRIT) 
DIMENSION  COEFt IDIM, IDIM) , XKAREA( IDIM) , RS ( IDIM) 
DIMENSION  fi ( 20 , ZO ) , XI ( 20 ) , Y ( 20 ) 

PRINT*," 

PRINT*," 

PRINT*," 

PRINT*,"  KXR  SOLN  FOR  ",N,"  TERM  EXPANSION  " 

C  PICK  METHOD 


I F 

<  I 

C.RIT 

.EG. 

,  1)  GO 

TO 

1 

GO 

TO 

20 

* 

X 

PR  I 

N»7 

*, " 

COLLOCATION  " 

CAL 

i 

CCU  COE? 

:  >  R  S  ,  N  , 

IDI 

K) 

20 

IF 

(i 

CRIT 

.  E3 . 

,2)  GO 

TO 

2 

GO 

TO 

30 

2 

PRI 

NT 

*  r  " 

GALL 

.ERKIN 

" 

CALL 

GAL < CGEr 

' , RS , N , 

IDIM) 

30 

IF 

(  I 

CR  -  7 

.  EG , 

,3)  GO 

TO 

3 

GO 

TO 

40 

3  PRINT*,"  LEAST  SQUARES  H 
CALL  LS ( CQEF ,R3 , N , IDIM) 

C  PRINT  RESULT  OF  CALL 


40  PRINT*,"  RIGHT  SIDE  VECTOR3  * 

PRINT*,"  CuEF  MATRIX3  " 

DO  5  J  3 1 , N 
DO  100  KK  =  1 ,  N 
100  A( J i KK ) 3  COEF(J.KK) 

Y ( J ) 3  RS(J) 

5  CONTINUE 

C  FIND  SOLN  CONSTANTS 
Ki-=1 
I D  G  T  3  8 

CALL  LEQT1F(C0EF,ML,N, IDIM.RS. I DGT, NX AREA, IER) 

PRINT*,"  IER3  " , IER 

IF  (IES.NE.34)  GO  TO  50 

PRINT*,"  EXACT  METHOD  SG-N  CONSTANTS3  " 

PRINT*, ("  A ( " , J , " ) 3  " ,RS< J) , J31 ,N) 

SORCON3  l.OE-B 
DC  GO  J  =  1 , N 
60  XI (J) 3  R S ( J  ) 
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-  1.0 

CALL  SCR  ( A  i  X I ,  RS ,  Y .  XKAREA , N , I D IN , ALP , SQRCON ) 

50  CONTINUE 
PRIM*  i  " 

PR  I  Ml  *  r SOLN  CONSTANTS  " 

PRINT*, ("  A  (  "  ,  J »  "  )  =  "ii?S(J),J  =  I,N) 

PRINT* , " 

RETURN 

END 

SUEREUT  I  \£  COL  ( CCEF ,  RS ,  N ,  I  DIM ) 

DIMENSION  COEF(  IDIM. IDIM)  ,RS( IDIM) 

C  FILL  RS  VECTOR 
DX=  1 . 0/ ( N+ 1 ) 

DO  1  j  =  1 .  N 

1  RS  <  J ) =  DX*  J 

C  FILL  CCEF  MATRIX 

DO  2  J -  1 1 N 
DO  2  I = 1 , N 
X  =  D  X  *  j 

CCEF ( J  .  I )  =  (I*(I+i )*(X**(I-1) > >+(X**(I+l) )-<X**I) 

1  -  <I*(I-l)*(X**(I-2))) 

2  CONTINUE 
RETURN 
END 

SUB ROUT  I  (\E  GAL  ( COEF ,  RS ,  N ,  I  DIM ) 

DIMENSION  CCEF ( I  DIM, ID IM) ,RS(IDIM> 

C  FILL  RS  VECTOR 
DO  1  J  =  1 » N 

1  RS(J)=  1 .0/ ( ( J+2 )  * ( J+3 ) ) 

C  FILL  CCEF  MATRIX 

DO  2  J  =  1  r  N 
DO  2  1=1 ,N 

COEF ( J , I ) =  (1.0*(((!+1)*I)+(I*(I-1)))/(J+I))  +  (2.0/1 J+I+2) ) 

1  -  (1.0/(J+I+3))-(1.0*I*(I-I)/(J+I-l))-(1.0*((I*(I+l))+l)/(J+I+l>) 

2  CONTINUE 
END 

SUBROUTINE  LS(CGEF,RS,N, IDIM) 

DIMENSION  COEF ( IDIM, IDIM), RS (IDIM) 

C  FILL  RS  VECTOR 
DO  1  J= 1 , N 

I  RS( J ) =  1.0  -  ( 1 . 0/ ( ( J+2 ) * ( J+3 )  ) ) 

C  FILL  CGEF  MATRIX 

DO  2  J  = 1 ,  N 
DO  2  1=1 fN 

IF  ((I+J).£Q.3)  GO  TO  11 
GO  TO  20 

II  COEF  ( J ,  I )  =  (1.0*((J*(J-1))  +  (J*(J  +  1)*I*(I  +  1))  +  (I*(I-1.')>/(1+J-1)) 

1  -  (1.0*( ( J*( J-l ) )  +  ( J*( J  +  l ! )  +  ! I*( 1  +  1 ) )  +  ( I*( 1-1 )  )  ) /( I+J)  ) 

2+  (1.0*<(I*(I+l))+i.0+(J*(J+l)))/(J+I+l))  +  ( 1 . 0/ ( I+J+3) ) 

3  -  (2.0/ ( I  +  J+2) ) 

4  -  (1.0*((J*U-l)*I*(I+l))+(I*(I-l)*J*(J+13))/(I+J-2)) 

GO  TO  SO 

20  IF  ( (I+JI.EQ.2)  GO  TO  21 
GO  TO  30 
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Zi  COEr  (  J  ,  I )  =  <:.0*(J*(J-!)*I*(I-!))/<I+J-3)>  -  (2. 0/ ( I+J+2 ) ) 

z  +  (i .o*( ( j*{ j-i ) )  m  j*o: )*i*( :+i ) )*<  i«(  i-i ) )  )/<  i+j-n  > 

3  -  ( 1 .0<-( !  J*(  J-  •  ) )  +  ( J*(„>: ) )  +  ( I*<  I*: )  )t(  :*(  2-1 ) ) )/( I* J? ) 

4  +  ( 1 . 0#  ( (  <  * ( I*  1 )  )-*-1.0+(j*(J  +  i) ) )  /  {  J  +  i  ♦  1 ) )  +■  ( 1 . 0/  ( I  +  j+3 ) ) 

GO  TG  90 

30  COEF  ( J 1 1 )  =  <1.0*(J*(J-l)*I»C-l))/<I  +  J-3))  -  ( 2 . 0/  ( I  +  J  +  2 ) ) 

2  -  (1 .  0*  (  (  J*  ( J-l  )  *  I*  (  I-rl )  )+(!*(  1-1 ) *J* ( J+l )  ) )/( I  +  J-2) ) 

2  v  {1.0*(<J*(J-1>)+(J*<J+1>*I*<I+1))+(I*<I-1))>/(I+J-1>) 

3  -  (1.0i((Ji(M)l  +  (J*(J  +  l))  +  (U(I  +  :))  +  (I*(I-l)n/(I+J)) 

-  -  (i.o«(c*(i+i))+i.o+(.j*(j+i)))/(j+i+i))  +  u.o/(!+j+3>) 

SO  CONTINUE 

2  CONTINUE 

RETURN 
END 

SUBROUTINE  COMPMRKRS.N, IDIM) 

D INENSIGN  RS ! I  DIM ) 

C  COMPARE  AT  MESH  POINTS/CCLLOCATIGN  POINTS 

PRINT*,"  ERROR  OF  MXR  AT  M E 3 n  POINTS  OF  FD  OR  COLLOCATION  POINTS" 
IT0L=0 
M  =  N 

DX=  1 .0/ (N’+i ) 

GO  TO  50 

30  IF  (ITOL.EQ.l)  GO  TO  31 

GO  TO  40 

31  M-99 

DX  =  .01 

PRINT*,"  ERRORS  FOR  WHOLE  REGION  AT  99  POINTS  " 

GO  TO  50 

40  IF  ( I TOL . EG . 2 )  RETURN 

50  ERRM=  0.0 
FERRM--  C.O 
XME=0.0 
XMF -  0.0 
DO  1  J  =  1  ,M 
X  =  DX* J 

Y=  ( S I N  ( X )  /  S I N  ( 1 . 0 ) )  -  X 
YM  =  0.0 
DO  10  I = 1 , N 

10  YM  =  YM  +  (RS< I )*( (X**I )-(X**(  1  +  1 ) ) ) ) 

ERR=Y-YM 

FERR=ERR/Y 

IF  (A3S(ERR) .GT.ERRM)  GO  TO  2 
GO  TG  5 

2  ERRN  =  ASS (ERR) 

XME  =  X 

5  IF  (A3S(FERR) . GT.FERRM)  GO  TO  3 

GO  TO  3 

3  FERRM=  ASS(FERR) 

XMF = X 

S  IF  ( ITOL.tQ.2)  GO  TO  1 

1  CONTINUE 

ITOL=  ITOL  +  ! 

PRINT*,"  MAX  ERI?OR=  “ ,  ERRM ,  "  AT  X=  " ,  XME 

PRINT*,"  MAX  FRACTIONAL  ERROR*  " , FERRM , *  AT  X=  ",XV- 
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50  TO  30 


SwO*. OL'TINE  -DSM C CGEF ,FD , HNAREA ,  N ,  I D I M ! 

DIMENSION  CGcFCD:K.:D:«),rD(IDIM>.kKA.?EA(IDIM) 

M=  K-l 

C  FILL  FD  VECTGR  AS  RIGHT  SIDE 
DX  =  1.0/N 
DO  1  J  =  1 ,  M 

1  FD<J)=  DX*J/ (N**2) 

C  FILL  CCEF  MATRIX 

DO  2  j  =  1 » M 
DC  2  1=1. M 

2  CCEF ( J . I ) =  0.0 

DIAG  =  2.0  -  <1.0/(N**2)> 

DO  3  J  =  1,K 

3  CGEF ( J  >  J )  =  DIAG 

IF  (N.EQ.2)  GO  TO  10 
M2=  N-2 
DO  4  J=2»M2 
KM= J-l 

COEF(JrKM)*  -1.0 
K.°=  J+l 

4  CCEF (  J  >  KP )  =  -1.0 
CCEF (1.2)=  -1.0 
CGEF ( M . M2 )  =  -1.0 

C  PRINT  RESULTS 
10  PRINT*," 

PRINT*," 

PRINT*," 

PRINT*,"  FINI'E  DIFFERENCE  SOLN  FOR  " » K , "  MESH  POINTS 
PRINT*."  RIGHT  SIDE  VECTOR*  ’ 

PRINT*, ("  RS ( " , j , " ) =  " ,FD( J) , J=1 ,M) 

PRINT*,"  CGEF  MATRIX*  " 

DO  12  J*!,M 

12  PRINT*, ( "  C(",J," ,1," )=  " ,COEF( J , I ) , 1=1 ,M) 

C  CALCULATE  SOLN  USING  THIS  IMPLICIT  SCHEME 

ML  =  1 
IDGT  =  8 

CALL  LEGT IF ( CCEF , ML , M> IDIM.rD, IDGT, NK AREA, IER) 
PRINT*," 

PRINT*,"  FINITE  DIFFERENCE  SOLN=  " 

PRINT*, ("  FD ( " , J , " ) =  " ,FD ( J ) , J= 1 , M) 

PRINT*,"  MESH  PGINT  LGCATICNS+  " 

DO  13  1  =  1,  M 

13  UK AREA ( I ) =  DX*I 

PRINT*,  ( "  :<(",I,")=  "  ,  WKAREA(  I ) ,  I  =  1  ,M) 

PRINT*," 

RETURN 

END 

SUBROUTINE  COM?FD(FD,N, IDIM) 

DIMENSION  FD(IDIM) 

C  COMPARE  ACCURACY  OF  FD  AT  MESH  POINTS 

PRINT*,"  ERROR  OF  FD  AT  MESH  POINTS=  " 

DX--  1.0/N 


•  • 
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M  =  N-l 

0.0 

FERRM 1=  0.0 
XME  =  0.0 
XMF*  0.0 
DO  1  J=1 »M 
X*  Da*  J 

y=  c  s : m  <  x ) / s : n ( i . o » )  -  x 

ERR*  Y-  FD(J) 

FERR*  ERR/ Y 

IF  (ABS(ERR)  .GT.ERRM)  GO  TO  2 
GG  TO  5 

2  ERRM*  ABS(ERR) 

XME*  X 

5  IF  (A3S(FERR) .GT.rERRM)  GO  TO  3 
GO  TO  6 

3  FERRIC*  ABS(FERR) 

XMF*  X 

G  PRINT*."  FOR  X*  *,X,"  ERROR*  ” .ERR,"  AND  FRACTION  ERROR*  " . FERR 
1  CGNTINUE 

PRINT*,”  NAX  ERROR*  "  , ERRM, "  AT  X*  ",X«E 

PRINT*,"  MAX  FRACTIONAL  ERROR*  " , FERRM , "  AT  X*  ",XMF 

RETURN 

END 

SU3R0UT INE  SGR ( A , XI , XO, Y, NX , N , IDIM, ALP , CONV) 

DIMENSION  A ( I D I M , I D I M ) , X I < I D I M > . Y ( I D I M ) ,WK(IDIM> 

i  ,xo( id in ) 

c 

C  SUCCESSIVE  OVER  RELAXATION  FOR  AX=Y 
C  XI  IS  AN  INITIAL  GUESS  AND  ALP  IS  A  SOR  PARAM 

C 

DO  1  I = 1 , N 
1  KK(I)*  XI(I) 

C  COUNT  ITERATIONS  FOR  MIN  GF  10  AND  MAX  OF  1000 
ITC*  1 

C  ALGGRYTHM  FOR  SOR 
B  DO  20  1*1, N 

XC  < I ) =  XX(I)  +  ( ALP*Y ( I ) /A ( I , I ) ) 

DO  9  J  =  I , N 

9  XO(I>=  XO(I)  -  (ALP*A(I,J)*WK(J)/A(I,I)) 

IF  (I.Ea.ll  GO  TO  20 

IM1*  1-1 
DO  10  J  —  1 , 1 M 1 

10  XO(I)=  XO(I)  -  ( ALr  *A( I , J )*X0( J ) /A( I , I ) ) 

20  CONTINUE 

C  CONVERGENCE  TEST 
I CONV*  0 
DO  50  I  =  1 ,  N 

IF  (A3S( ( XO ( I )-WK ( I ) )/X0( I ) ) .GT.CONV)  ICONV*  1 

wxm*  xom 

50  CONTINUE 

I-  ( ICCNV.E3.0)  GO  TO  80 
IF  ( ITC. GE. 1000)  GO  TO  70 
ITC*  ITC+1 
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sc;?  NOT  CONVERGED 


GO  TO  8 

70 

RETURN 

80  IF  C7C.GT.iO)  GO  70  90 
ITC=  I7C+1 
GO  TO  S 

SO  PRINT* ,  “  SCR  CONVERGED  * 
RETURN 
END 
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Transient  Computer  Program 


PRCCRA^  TPEi  ( INPUT  i  OUTPUT » ) 

CC1P-EX  AL“A ( 20 ) ,2(20,20) .LAMDA(ZO) 

REAL  A  ( 20  >  ZO )  ,3(20,20)  AS  (800)  ,  RZ<  £00 ) ,  RALFA  ( 40 ) ,  BETA  ( 20) 

1  , RLAMDh(ZO)  ,.332(20,20)  ,FD(  100, 2) 

ECU!  VALENCE  <  AlFA  ( 1 )  ,  PALP  A  ( 1 ) )  ,  ( Z  ( 1 , 1 ) ,  RZ<  1 ) ) 

1011=  20 
10112=  2*IDIM 
IDI.V.S  =  IDI>!*IDIM2 
ID  1 MF=  100 
DO  5  j  J  =  1 , 3 
T  =  .05  *  JJ 
DO  !  ICRIT=1,2 
DO  1  N  =  1 » 20 

CALL  MwRS0_( A,S ,N»  I  DIM,  I CR IT, Z , RRZ.RLAKDfi , ALFA , BETA , KK ,RZ , 
1  RALFA, IDI12, ID; IS ,  LAND A) 

CALL  KWRCCMCRRZmRLAKDA,  IDIK,N,T) 

1  CONTINUE 

DO  3  11=1,10 
N=  11*10 
DX=  1 . 0/ (N+l ) 

DT  =  (DX**2)/2. 1 
IF  (DT.CE.T)  GO  TO  10 
IVAL=  T/DT 
I VAL=  IVAL+1 
DT  =  T/IVAL 
GO  TO  30 
10  DT=  T 

C  LOAD  INITIAL  VALUES 
30  DO  4  1  =  1,  N 

4  FD ( 1 , 1 ) =  0 . 0 

PRINT*," 

PRINT*," 

PR  IN’*," 

PRINT*,"  FD  SOLN  FOR  ",N,"  MESH  POINTS  n 
DO  2  J  =  1 , 1 VAL 

CALL  FDSGLN(FD,AA,N,DT ,HX, IDIMF) 

2  continue 

PRINT*, CM  FD ( " , I , " ) =  " ,  FD ( 1 , 2 ) , I  =  1 , N) 

CALL  CCMPFD(FD,N, ID  IMF, T) 

3  CONTINUE 

5  CONTINUE 
END 

SUBROUTINE  GALT ( A,3,N> IDIM) 

DIMENSION  A ( ID IDIM)  ,8 (IDIM,  IDIM) 

C  FILL  A  AND  3  MATRICES 
DO  1  J=1,N 
DO  1  I = 1 , N 

.  A(J , I ) =  (<1.0#(I*(I-1)>)/(J+I-1>) 

1  -  ((2.0*(I*+Z))/(J+I)) 

2  +  < ( 1.0*(I*( 1  +  1 ) ) )/( J  +  I  +  l ) ) 

B ( J  ,  I  >  =  (2.0/C J  +  I+Z) )  -  ( 1 . 0/ ( J+I  +  l ) ) 

•  • 
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1  -  ( 1 .0/( J  +  I+3) ) 

i  CONTINUE 
RETURN 
END 

SUBROUTINE  COLT ( A , B , N , I  DIM ) 

1)1'  ENSIGN  A<  IDIM,  IDIM)  ,3<  IDIM,  IDIM) 

C  FIlL  A  AND  B  MATRIX 
DX-  1 .0/ (N+ 1 ) 

DC  i  J=1.N 
DO  1  1  =  1,  N 
',{■■■  DX*  J 

I-  (I.EO.l)  GG  TO  10 
GO  TO  20 

10  A  ( J ,  I )  =  -2.0 
GO  TO  50 

20  IF  (I.EQ.2)  GO  TO  30 

GO  TO  40 

30  A <  J  > I ) =  2.0-(6.0*X) 

GO  TO  50 

40  A ( J . I )  =  <I*(I-i)*(X**(I-2)))  -  ( I*(  1  +  1 )*(X**(  1-1  )  > ) 

50  S ( J , I ) =  (X**(I  +  1>)  -  ( X  *  *  I ) 

1  CONTINUE 

RETURN 
END 

SUBROUTINE  MNRSCL(A,3,N. IDIM, ICRIT,Z,RRZ,RLA«DA, ALFA, BETA, 

1  XK,R?.,RALFA,  I D I M2 , 1 D 1 MS  >  LAMD A ) 

COMPLEX  AL.rA(  IDIM)  ,Z(  IDIM,  IDIM)  ,LA.V.DA(  IDIM) 

REAL  A ( IDIM, IDIM) , 3( IDIM, IDIM  >, WK ( IDIM5 ) ,RZ( IDIMS) , RALFA ( IDIM2 ) , 
1  8E7A(  IDIM)  ,RLAi'DA(  IDIM)  ,RRZ<  IDIM,  IDIM) 

C  COMPUTE  A  AND  B  MATRICES 
PRINT*," 

PRINT*," 

PRINT*," 

PRINT*,"  MNR  SOLUTION  FOR  ",N."  TERM  EXPANSION  " 

IF  ( ICR  IT. EG.  1 )  GO  TO  10 
GO  TO  20 

10  PRINT* , "  COLLOCATION" 

CALL  COLT <A,B,N, IDIM) 

GO  TO  40 

20  IF  ( ICR  I T .  EG .  2 )  GO  TO  30 
GO  TO  40 

30  PRINT*,"  GALLERKIN" 

CALL  GALT(A,B,N, IDIM) 

40  PRINT*,"  MATRIX  A=  " 

DO  1  J  =  1 >  N 

1  PRINT* , ( "  A (",J,", ",!,")=  " , A( J , I) , 1=1 , N) 

PRINT*,"  MATRIX  B=  " 

DO  2  J  =  1 ,  N 

2  PRINT*, ( ”  B(",J,",",I,")=  ",B(J,I),I=1,N> 

C  CALCULATE  EIGENVALUES  AND  EIGENVECTORS 

I JOB=  2 

CALL  EIGZFIA, ID IM , B , IDIM.N, I JG3 , RALFA , BETA , RZ , I  DIM, NX, ICR) 
PRINT*,* 

PRINT*,"  SOLUTION  VALUES  " 
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print*,"  i Er?= 

3'.:  s '■>, performance  index3  mmi) 

C  .-U’M.  (.-..f"i_EX  AND  EcAl  cIGENVALutS 

L’U  w  J  -  1  ,  N 

I."  ( EETAf  j  ) .  EG .  0 . 0 )  GO  TO  100 
0G  TO  i 10 

100  PR I  NT*  ? "  SETA ( " . J . “ ) =  0.0  " 

ETCP 

110  LA>DA(J)=  ALFA< J>/3ETA< J) 

R-AMDA ( J ) 3  ALTAI J ) /3ETA( J ) 

3  CONTINUE 

PRINT*,"  COMPLEX  EIGENVALUES3  “ 

PRINT*, ( "  _C(U,I,“)-  " , LAMDA ( I )  •  1 3 1 ,  N ) 

PRINT* REAL  EIGENVALUES3  " 

PRINT*, ("  LR(",I,")=  “ ,RLAMDA( I ) , 1= 1 ,N) 

PRINT*,"  COMPLEX  EIGENVECTORS3  " 

DO  4  J  =  J ,  N 

4  PRINT* ,  ( "  EV ( " >  J , " , " , I , " )  3  ",Z(J,I),I  =  1,N) 

C  FORM  REAL  EIGENVECTORS 

C AoL  CXXF IX ( Z , RRZ , N , I D I M , RZ , I D I  MS ) 

C  SETUP  MATRICES  TO  CALCULATE  EIGENVECTOR  MULTIPLIERS 
IF  .( ICRIT.E3.1)  GO  TO  200 
GO  TO  210 

200  CALL  C3LV(A,3ETA,RRZ,N, I  DIM) 

GO  TO  230 

210  IF  CCRIT.E3.2)  GO  TO  220 
GO  TO  230 

220  GALL  GALVt A , BETA , RRZ ,N , IDIM) 

C  SOLVE  SYSTEM  FGR  MULTIPLIERS 
230  ML  =  1 

IDGT =8 

CALL  LEGT 1 F  ( A ,  ML ,  N ,  I  DIM! ,  SETA  ,  I DGT ,  KK  ,  IER ) 

PRINT*,"  FGR  I.C.  SOLUTION  IER3  ",IER 
PRINT*,"  MULTIPLIERS  ARE3  " 

PRINT*, ("  M ( " , I , " ) 3  “ , SETA ( I ) , I  3 1 >  N ) 

C  SCALE  EIGENVECTORS  BY  MULTIPLIERS 
DO  G  K 3 1 ,  N 
DO  G  I  3 1 , N 

B  RRZ ( I , K ) =  RR2( I , K ) *SETA (K ) 

RETURN 

END 

SUBROUTINE  GALV ( A , BETA , RRZ , N , IDI M) 

DIMENSION  A( IDIM, I  DIM) .BETA (IDIM) ,RRZ( IDIM, IDIM) 

C  FILL  BETA  MATRIX 

DC  1  J  =  1 ,  N 

1  B E T A ( J ) 3  (2.0/( J+2) )-( 1 .0/< J+l ) )-( 1 .0/( J+3) ) 

C  FILL  A  MATRIX 

DO  2  J31,N 

DO  2  K=1,N 

A ( J , K ) 3  0.0 
DO  2  1 3  1 ,  M 

A ( J , X ) 3  A( J,X>  +  (RRZ;Z.K>*((1.0/<J  +  I  +  l>>-(2.0/<JM+Z)) 
1  + ( 1 .0/ ( J*I+3) ) >  > 

2  CGNTiNUE 
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RETL'R.V 

END 

SUBROUTINE  CCLV  ( A .  BETA ,  RRZ , N,  I  Dll". ) 

DIMENSION  A(  IDIM, I D 1“ ) , BETA ( I  DIM ) ,RRZ( IDIM, IDIM) 

dx=  i.o/<n+i> 

C  FILL.  SETA  MATRIX 
Du  1  J  = 1 >  N 

1  £["A( J ) =  ( D  X  *  J  )  - 1 . 0 
C  FILL  A  MATRIX 

DC  2  J - 1 ,  N 
X -  DX*J 
DO  2  K  =  1,N 
A< J  r  K  )  =  0.0 
BO  2  I  = 

A(  J  r  K )  =A  (  J,K)+(RRZ(  I  ,K)*(  (X**I  )-(»*(  1  +  1) )  > ) 

2  CONTINUE 
RETURN 
END 

SU3RCUT I N'E  SOLNT  ( X ,  T ,  TS ,  CONV ) 

PIE  =  3 . 1 4 15S2S53G 
TSG=  1.0-X 
11=1 

1  E>!=  (?1**2)  »(PIE**2)*T 

ARG  =  X*PIE*X 

TS"  TSO  -  (2.0*3IN(ARG)/(?IE*.M*EXP(EX) ) ) 

IF  ((2.0/(P:E*K*EX?<EX))).L7.CO.W)  RETURN 

3  TSO=  TS 
M  =  (1+1 
GO  TO  1 
END 

SUBROUTINE  *XRST(X,7,T(1,RRZ,RLAMDA,  IDIM,N) 
DIMENSION  RRZ( IDIM. IDIK) ,RLANDA( IDIM) 

TM=  l.O-X 
DO  1  1=1, N 
CT  =  O.O 
DO  2  J  =  1  f  N 

IF  ( RUANDA < J) .IT. O.O)  GO  TO  2 

IF  ( (RLANDA(J)^T) .GT. 700.0)  GO  TO  2 

CT  =  CT+(RRZ(IrJ) /EX? ( R _Ai1DA ( J ) *T ) ) 

2  CONTINUE 

7K=  TX+(C7<M (X**I >-<X**< 1+1) ) )  ) 

1  CONTINUE 

RETURN 
END 

SUE  ROUT  I NE  CVRCOM  (  RRZ .  RUANDA ,  I D I M ,  N ,  T  ) 

DIMENSION  R.RZ(  IDIK,  IDIM)  .ELAMDAI  IDIM) 

CONV=  1.0E-20 
C  COMPARE  AT  -D  POINTS 
DX  =  1 . 0/ (N+ 1 ) 

ERRM=  0.0 
FERRK=0.0 
XMF-0.0 
XME=0.0 
DO  1  1=1, N 
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a  ■  D  X  *  I 

r  ,LL  s:^\to:,t,te,conv> 

L.  —  A  T-  '  ■>  i  ,  i  "  r  it  R ,  .r  _  A.'  -J  A  ,  i  i‘  -  I  N  ) 


FEaa=  E  \’  R  /  T  S 

IF  <A3S(ERR)  .G7.ERRM)  GO  TO  10 
GO  TO  TO 

10  E  R  R  M  -  n3S  (  ERR ) 

XMc  -  X 

20  IF  ( A30  ( FERR  ) .  GT .  CERRM )  GO  TO  30 
GO  TO  1 

30  FERRM*  AES  (FOR?) 

xmf=  x 

1  CONTINUE 
PRINT*," 

PRINT*,"  ERROR  FOR  FD  POINTS  AT  T=  ",T 
PRINT*,"  FAX  ERROR1  " , ERRM , "  AT  X*  ",XME 
PRINT* , 11  Ffi;.;  -RACTIuNA.  ERROR-  " ,  FERSF ,  “  AT  X*  ",XKF 
0  COMPARE  AT  S3  POINTS 
DX=  .01 
ERRi1  =  0 . 0 
FEiRRM-0 . 0 
XNF-0 . 0 
XNE-0.0 
DO  2  1=1,99 
X=  DX*I 

CALL  SCLNT (X > T, TS , CCNV ) 

GALL  MARS  7 ( X , 7 , 7M , RRZ , RLAM3A , IDIM, M ) 

ERR*  TM-T3 
FERR*  ERR/TS 

IF  (ASS (ERR) . GT.ERRM)  GO  TO  50 
GO  TO  GO 

50  ERRM*  ABS(ERR) 

XMIE*  X 

GO  IF  (A33(F ERR) .  G7.FERRM)  GO  TO  70 
GO  TO  2 

70  FERRM*  ASS (FERR) 

XMF*  X 

2  CONTINUE 
PRINT*," 

PRINT*,"  ERROR  FOR  S9  POINTS  AT  T=  ",T 

PRINT* , "  MAX  ERROR*  " , ERRM , "  AT  X=  ",XME 

PRINT*,"  MAX  FRACTIONAL  ERROR =  ", FERRM, "  AT  X*  ",XMF 

RETURN 

END 

SU3RCUTI KE  CKN-  IX  ( 2 ,  RRZ ,  N ,  IDIM ,  R2 , 1  DIMS ) 

COMPLEX  Z( IDIM, IDIM) , CXI 
REAL  RRZ ( IDIM, IDIM) ,RZ( ID  IMS) 

C 

C  STORE  EIGENVECTORS  AS  REALS  AND  C:-'ECK  FOR  COMPLEX 
C  EIGENVECTORS  (XHICH  COME  IN  CONJUGATE  PAIRS)  AND 

C  SET  ONE  EIGENVECTOR  TO  THE  REA_  PART  AND  THE  OTHER 

C  TO  THE  IMAGINARY  PART  (NOIE  COMPLEX  EIGENVALUES  ARE 
C  TAKEN  AS  THE  ‘REAL  PART  FOR  APPROXIMATION ) 
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L. 

:s;<ip=  o 

DC  :  jm.\ 

IF  ( ISKIP.E3. 1 )  GO  TO  4 

:  ck  =  o 
:c  2  :=i,m 

RRZC.J)  =  Z(I.J) 

cti  =  zc.j)  -  rrzii.j) 
c  c:~ec\  for  complex  eigenvectors 

IF  (OKI .NE. (O.0. 0.0) )  ICK  =  1 

2  CENT  I. NOE 

IF  (ICN.EG.O)  GO  TO  1 
C  KATE  IMAGINARY  PART  NEXT  EIGENVECTOR 
C  AND  SKI?  TO  FOLLOWING  EIGENVECTOR 
DO  3  1=1 rN 

I  OFF  =  2  +  <  C-1  +  CDIK*<  J-l)  ))*2) 

RRZ ( I , J+ 1 ) =  RZ(IOFF) 

3  CCNTiNuE 

:skip=  i 

GO  TO  1 

4  I S  K I F  =  0 

1  CONTINUE 

PRINT*,"  REAL  EIGEVVECT ORS=  " 

DG  10  I  =  1  > N 

10  PRINT* , ( "  EV(",I,H,",J,")=  " » RRZ ( I »  J ) , J  =  1 , N ) 

RETURN 
END 

SUBROUTINE  COPPED (FD.N, IDIXFrT) 

DINENSIC:  FDC  IDIP1F.2) 

C  CCPFA'-E  SOLUTION  AT  PESH  POINTS 

D!<-  1 .0/ ( N+ 1 ) 

CCNV=  1. 02-10 
ERR,P  =  0.0 
F  E  R  R  X  =  0.0 
X,yZ=  0.0 
XNF=  0.0 
PRINT*," 

PRINT*,"  FD  ERROR  COXPARISION  AT  MESH  POINTS  " 
DO  1  I  =  1  ,  N 
X=  DX»I 

CALL  SG-NT ( X , T , TS , CONV > 

ERR=  FD ( I , 2 )  -  TS 
FERR=  ERR/TS 

IF  ( A3S ( ERR ) . GT . ERRM )  GO  TO  10 
GO  TO  20 

10  EERP= ASS ( ERR ) 

SCLPEF=  AES(FERR) 

XiyE  =  X 

20  IF ( A3S ( FERR ) .GT.FERRX)  GO  TO  30 

DO  TO  1 

30  FERRP=  AFS(rERR) 

SOLKFE=  ASS (ERR) 

X  X  F  =  X 

I  C  J  N  •  I  i\  'J  t. 
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